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20.  The  fpllowirg  transient  phano— ns  occur  In  a few  milliseconds : 

penetration  of  hot  8s***  Into  interstitial  voids,  convecciva  heating  of 
pallets  to  ignition,  granular  bad  compaction  and  rapid  pressurization  in 
a thlck-valled,  steal  chamber.  Pressure  transients  and  flaae-front  speed* 
are  measured  by  high-frequency  pressure  transducers  and  Ionization  probes, 
respectively.  A gaseous  pyrogen  ignition  system  utilizing  a spark  ignition 
was  found  nose  suitable  for  varying  igniter  strengths  and  achieving  high 
reproducibility. 

> 

The  physlo-cheelcal  phenomena  described  above  have  been  formulated 
into  a theoretical  model.  The  method  esq>loyed  in  this  stidy  is  the  gas 
dynamic  approach  which  is  developed  by  formulating  the  governing  equations 
on.  the  basis  chat  mass,  women  turn  and  energy  fluxes  are  balanced  over 
control  volumes  occupied  aeparataly  by  the  gas  and  partlcle-phasaa . The 
governing  equations  were  derived  in  the  form  of  coupled,  non-linear,  in- 
homogeneous partial  differential  equations.  This  system  of  equations  Is 
solved  together  with  soma  empirical  input  such  as  heat  transfer,  flow 
resistance.  Intergranular  stress  and  the  burning  rate  law.  They  are 
found  to  be  of  hyperbolic  type,  since  all  the  eigenvalues  are  real. 

The  compatibility  relations  are  used  as  extraneous  boundary  conditions 
together  with  some  ptayrJ  \1  boundary  conditions  to  solve  the  boundary 
points.  A stable,  fast  - onve rg an t , explicit  numerical  scheme  incorporated 
with  predictor  and  corrector  are  used  in  the  interior,  and  the  method  of 
"-^characteristics  is  used  on  the  boundaries. 

The  results  show  that  the  flame  front  eccnleretee  significantly  end 
the  reta  or  pressurization  increases  substantially  in  the  downstream 
direction.  The  igniter  strength  and  tba  propellant  gasification  temperature 
ware  found  to  have  pronounced  affects  on  the  pressurization  process.  The 
theoretical  predictions  are  In  close  agreement  with  the  experimental  data. 

rl 
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nmoDOCTioK 

Tba  activation  and  the  overall  objective  for  the  present  invencigatioa 
is  first  defined.  Having  provided  the  motivation  and  objective  for  this, 
aoaa  useful  applications  are  given  is  regard  to  this  study.  We  then  provide 
sues  background  information  x elated  to  pioneering  work  in  the  early  stage 
of  interior  ballistic  modeling.  Inference  is  ends  to  ease  of  the  related 
v*ork  in  this  area.  Ibis  chapter  is  concluded  by  e brief  discussion  of 
the  basic  approach  taken  for  different  models  associated  with  flame 
spreading  and  granular  propellent  bamirg.  the  mala  purpose  of  this  die-- 
cum  ion  la  to  point  out  the  fundamental  differences  in  the  method  of 
approach  for  all  noeals. 

x.l  Koclvsrt*.  and  Objective 

u 

'this  research  is  a dan  1st  ad  by  the  potential  of  gas— permeable/  pro- 
pellents for  achieving  high  gasification  rates  «d  producing  high  thrust* 
during  aactramaly  short  time  intervals.  The  research  ia  also  partially 
motivated  by  the  recant  progress  is  the  area  of  granular  propellent  studies 
(.1,2,3,*,  5).  Besides  the  encouragement  offered  by  these  advance- 
ments la  this  field,  the  problem  of  flams  spreading  is  porous  propellent 
chargee  ia  obviously  important  la  the  design  aad  analysis  of  propulsion 
systems. 

The  overall  mstrek  motivation  aed  objectives  ere: 

1)  Aeameloua  burning  behavior  hes  been  observed  ia  cocsPesdsn  of 
ttaakr  solid  ^repellents.  k batter  understanding  of  the 

*Tha  word  "gee-permeable  prey  11  ant**  or  porous  propellant"  used  In 
this  s'mdy  deal  nates  either  e pacfltnd  bed  of  granular  pellets  with  or 
withwu  perforation,  or  e cant  monolithic  propellant  with  perse. 

t 


transient  cuabustlon  phenonsna  la  required  Co  raduca  haaarda 
ard  to  laprove  the  design  of  propulsion  aystuss. 

2)  To  advance  tha  «tat*— of-cha-art  in  tha  conbuation  of  granular 
propellants  by  formulating  a coaplate  thaoraelcal  nodal  dascribJ.ufl 
tha  Important  phyaieal  phanonanw  involved  in  both  gaseous  and 
•olid  phaaaa. 

3)  To  solve  tha  theore  leal  nodal  by  selecting  a a tab la,  faat 
convergent  numerical  achana  ao  than  affacta  of  a)  Igniter 
strength  b)  propellant  phyaio-chani cal  charact aria tic  and 
c)  propellant  loading  density  can  ba  aeudiad. 

4)  To  conduct  aeparlnentaJ  firing*  in  a cfcick-ualled  eoabuator 
with  r.  fixed  boundary,  to  naaaura  tha  preaaure-tina  transients 
at  various  axial  locations  and  flana  propagation  ratac  ao  that 
tha  thaoraelcal  nodal  developed  in  this  study  can  ba  verified. 

Tha  results  obtained  fron  this  study  are  directly  uaaful  for  nore 
accurate  predictions  of  tha  interior  ballistics  of  artillery  ays tana.  Tha 
research  it  also  pertinent  to  the  futvra  duvelopaamt  of  and-bu  ruing 
porous  propellant  charges  to  be  used  in  a nora  prciictabla  and  con- 
t Tollable  for  rocket  or  ala  alia  propulsion.  In  addition,  the 

techniques  developed  in  this  investigation  can  lead  toward  both  a batter 
understanding  of  porous  propellant  combustion  aud  the  aarller  behavior 
of  the  overall  deflagratlon-to-detoaatlon  transition  (DDT) . 

The  pbaacana  that  occur  daring  tha  trmlant  arrant  include  tha 
folloarlng: 

a)  the  frsaatretlon  of  bet  p jduct  gases  into  tha  grannlar  bad 

b)  convective  beating  of  propellent  to  Ignition 

c)  effect! an  of  the  greasier  bed 


Id 


d)  generation  of  peak  pressure  vithin  the  chaaber 

e)  propegetion  of  the  peek  pressure  to  the  sheer  disc  end 
of  the  granular  bed 

f)  flaae  front  acceleration 

g)  rupture  of  shear  disc  which  initially  seals  the  combustion 

h)  the  depressurization  processes  after  the  burst  of  tho 
shear  disc 
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1.?.  Pioneering  Work 


The  theoretical  formulation  of  granular  propellent  combustion  for 
interior  ballistic  calculations  can  be  categorised  into  three  separate 
stages:  1)  lumped  parameter  analysis;  2)  theoretical  modeling  for 

earlier  behavior  of  the  ignition  and  flame  spreading  processes  by  assuming 
that  the  granular  propellents  are  fixed  in  space;  and  3)  modeling  cf 
mobile  granular  beds  so  that  the  overall  transient  processes  in  ballistic 
cycle  study  can  ba  achieved. 

Baer  and  Frankie  (6) were  the  pioneers  to  use  the  lumped  parameter 
method  to  develop  an  interior  ballistic  model.  Their  model  is  governed 
by  five  simplified  equations.  Due  to  the  uncertainties  of  the  burning 
rate  law,  this  model's  predictive  capability  is  limited.  This  is  further 
improved  by  introducing  ballistic  fitting  into  the  burning  rata.  Soma 
other  lumped  parameter  analyses  are  contained  in  Corner's (7)  text  including 
the  classical  work  of  Lagrange,  Love,Pidduck  *nd  Kant  on  hydrodynamic 
theories  of  interior  ballistics. 

It  was  not  until  recently  that  K.uo,  Vichnavetsky  & Summerfiald 
(l,  8) have  developed  a fixed-bed  theoretical  model  (referred  as  KVS  me .el 
in  ballistics , to  describe  flame  propagation  in  a packed  bed  of  granular 
propellent*.  Tha  calculated  pressure-time-distance  transients,  wave 
propagation  spaed,  and  mass  fraction  of  propellant  burned  during  flame 
spreading,  ell  agree  wall  with  tha  experimental  data  obtained  by  Squire 
end  Devine  (9)  in  e venting  bomb  experiment,  Soper(io)  conducted 
experimental  firings  in  soma  artillery  systems.  Ha  utilized  translucent 
cartridge  cases  together  with  flash  X-ruy  and  photographic  tachnlques  to 
detazalne  tha  relative  positions  of  tha  ignition  front  to  tha  pressure 
wave  front  in  the  packed  bad  of  NACO  propellents.  Be  claims  that  hi' 
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experimental  data  also  agrees  wall  with  the  theoretical  conclusion  made  by 
chs  KVS  modal  (j  , 8)  i.e.,  tbs  Ignition  front  liss  near  tbs  leading  sdge 

of  tbs  prsssurs  wave.  Kuo  st  al  (l,  8)  used  an  implicit  finlts  difference 
schema  to  obtain  tbs  numerical  solutions  and  wars  also  checked  by  Kitchens 
(11) using  chs  method  of  characteristics. 

The  TVS  model  and  its  numerical  solutions  (1,  8)  basud  upon  the  fixed 
bed  assumption  for  a packed  bad  of  granular  propellents,  has  stimulated 
the  lntarast  of  other  investigators  in  pursuing  studies  in  both  theoretical 
and  experimental  directions.  Governing  equations  similar  to  the  KVS  model 
equations  have  been  further  developed  by  several  investigators,  to  include 
the  motion  of  solid  propellant  grains.  The  model*  ere  developed  by  Fisher  an* 
Graves  (2),  Gough  and  Zwerts(3;,  East  and  McClure  (4),  Krler,  Van  Tassall, 

Raj  an  end  VerShaw  (5jt  Gough  (l?),  and  Kuo  and  Summer  field  (30(Scuw  discussions 
on  each  modal  are  shawm  in  sac cion  2.1). 

7a  the  formulation  of  boundary  conditions  for  e transient  combustion 
study,  axcapt  KVS  model  (8)snd  Gough's  work  (l^),  none  of  the  previous  in- 
vestigators have  used  the  compatibility  relations  for  the  solution  of 
boundary  conditions . The  treatment  of  boundary  conditions  are  usually  not 
satisfactory.  Dus  to  tha  complication  «f  the  theoretical  modeling  and 
difficulties  In  experimentation,  no  investigators  have  verified  their  models 
by  performing  their  own  experiment*  to  match  with  the  same  conditions  for 
tbs  models. 

la  the  present  investigation,  experimental  firings  axe  done  to  verify 
our  theoretical  model.  The  governing  equations  are  derived  using  the 
pr  ftciples  of  gee  dynamics.  The  nature  of  the  system  of  governing  equations 
ere  determined  to  be  of  hyperbolic  type.  The  compatibility  rolitfiona  are 
deduced  to  supplement  the  physical  Boundary  conditions  to  solve  the  fe«  ucidety 
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point*.  A combination  of  explicit  finite  difference  scheme  atm  the  method 
of  characteristics  is  sstployad  for  the  numerical  solutions.  The  pre- 
dictions made  by  the  present  model  are  found  to  have  good  agreement  with 
our  experimental  date'. 


Dut  to  Che  complexity  in  describing  the  physical  processes  for  chi* 
two-phase  combustion  problem,  •ever*',  different  approaches  have  been  pro- 
posed. The  theoretical  ate chocs  c^a  b*  classified  into  four  categories: 

A.  statistical  rod el ; 

B.  continue  mechanics  models; 

C.  formal  averaging  models;  and 

D.  two-phase  fluid  dynamic  models 

The  method  employed  la  this  study  is  the  two-phase  fluid  dynamic  model 
which  is  developed  by  formulating  the  goweraix-^  equations  cn  thj;  basis 
that  mass,  momentum  end  energy  fluxes  are  balanced  over  control  volumes 
occupied  spearataly  by  the  gee-  and  particle-phases. 

Dp  to  now,  statistical  methods  have  not  been  applied  to  this  particular 
combustion  problem  due  to  their  Inconvenience  in  obtaining  useful  resvlts 
end  also  due  to  the  luck  of  reliable  statistical  information  for  describing 
some  of  the  important  terms  in  the  model.  Only  a few  applications  of 
these  methods  are  found  with  application  to  the  governing  equations  cf 
chemically  non-reacting  tvo-phesc  flow  system.  statistical  methods 

ere  included  in  the  text  by  Be  ran  (13),  A'ong  all  the  investigators , such 
as  Aronov  (j.*^  Law  and  Fung  ^.5),  Cox  and  Brenner  the  moat  comprehensive 
vorks  are  those  dona  'rj  Buyer ich(l7, 18 , 19).  The  main  difficulty  of  this 
approach  lias  in  the  description  and  definition  of  the  mean  properties  of 
e twpbcsi  motion  and  also  in  the  evaluation  of  the  net  effect  of  micro- 
cr.opic  behavior.  Because  of  these  difficulties  involved,  the  statistical 
approach  is  not  followed  in  the  present  investigation. 

In  the  continuum  mechanics  methods,  the  gsa  and  aolid  media  are  con- 
sidered at  separata  can  inua  that  interact  with  each  other  based  on  the 


continuum  mechanics  principle*.  Tha  basic  theory  was  laid  down  by 
Truewdall  and  Toupln  (20),  Most  analysis  In  these  methods  start  with  the 
conservation  lavs.  Since  continuum  mechanics  dous  not  seek  to  reconcile, 
the  properties  of  matter  with  the  structure  of  the  constituents,  therefore, 
this  theory  has  been  applied  both  for  hoocgeneoua  fluids  and  heterogeneous 
mixtures , 

Earlier  studies  by  Green  and  yagh-ai  (21,  22),  Erigan  and  Ingram,  (23)  and 
othars  ware  criticized  by  Trueadall  (24)  .Muller  (25)  and  Porl*(26)  have 
investigated  the  constitutive  lavs  for  mixture  of  tvo  fluids.  Hone  of  thu 
previous  invasti, gators  have  made  jn  attempt  to  conatruct  a theory  with 
immediate  application  to  the  tvo-plmse  reacting  flow  coubuction  problem. 
Recently,  Krier  et  el (5)  based  on  the  continuum  uechenics  concepts,  formu- 
lated equation*  for  a modal  that  associated  with  flame  spreading  aud  com- 
bustion in  porous  propellent  beds.  Ths  weakness  of  the  model  developed  by 
Krier  at  al  la  discussed  In  section  2.1. 

The  formal  averaging  methods  ausume  Chet  the  f J ow  behavior  of  the  micro- 
scopic aspects  JLs  detorminrd  snd  cask  for  the  macroscopic  aspects  of  the 
flow,  the  formulations  are  derived  by  averaging  the  flow  properties  over 
regions  large  compered  with  the  size  of  propellant  grains,  but  small 
when  coopered  with  the  macroscopic  scale  of  the  experiments . Among  thou* 
investigator®  who  applied  these  methods,  the  bectsr  known  works  are  three 
by  Slattery  (27,  2%  Anderson  hid  Jacsaon  (2$jt  Whitaker  (3(J,  and  Fauton  (7V 
Recently.,  Gough  (1%  based  on  the  works  by  Arderaon  and  Jadkaon  (29i  Wbiuakar 
(H&  a ru  P ancon  (33)  developed  a theoretical  model  r.o  solve  the  fiw  of  t 
coMpresaible  gas  through  an  aggregate  of  mobile  reacting  perticlai . The 
differences  in  derivation  of  conservation  aquations  between  our  approach 
and  tha  formal  averaging  approaches  era  analogous  to  tha  differences 
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bstvsan  the  Boltrmaim'c  Molecular  derivation  (32)  and  Cha  conventional  con- 
tinuum derivation  of  conservation  aquations  In  fluid  aachanics . Tha 
difference  batwaen  thaaa  two  approach**  are  discussad  in  mors  detail  in 
Chapter  11. 

Filially,  tha  evo-phase  fluid  dynamic  methods  are  considered  by 
formulating  tha  governing  equations  on  tbs  basis  that  mass,  momentum 
and  energy  fluxes  ere  balanced  over  control  volumes  occupied  separately 
by  the  gas-  and  particle-phases.  Among  the  previous  investigators,  Murray 
(33),  Kraiko  and  Sterain  (34),  Kigma-.ulin  (35),  East  and  McClura(4)  , Fisher  and 
Craves  ( 2),  and  Kuo  and  iiinmserfiold  (3$  have  formulated  complete  vets  of 
governing  equations  for  the  two— phase  reacting  flow  problem.  The  Interface 
conditions  between  the  gas  and  the  solid  are  considered  by  local  continuity 
of  fluxes  using  empirical  correlations . Tha  advantage  of  using  cjiis 
approach  is  to  be  able  to  lump  the  detailed  flow  behavior  between  tha  wo 
itUnartu  into  tarns  which  ren  b*  evaluated  by  using  empirical  correlations 
such  es  heat  transfer,  flow  resistance.  Intergranular  stress  and  burning 
rate  Lav. 

Mathematically,  this  problem  la  .nporoachad  by  deriving  tha  mass, 
mount  im  and  energy  equations  for  tbs  gas  phase  and  the  mass  and  momentum 
equations  for  the  solid  phase  in  s transient  quasl-ooe-dimecaional  form. 

In  summary,  statistical  methods  are  rejected  because  in  this  investi- 
gation ve  ere  interested  in  the  mesa  values  end  there  is  not  enough  in- 
formation to  evaluate  the  net  contribution  of  the  microscopic  fluctuation 
term.  The  continuum  mechanics  methods  era  not  considered  because  the  co- 
existent* of  solid  end  gee  are  resumed  without  considering  tbs  inter- 
penetrating between  the  two  media.  The  formal  averaging  methods  ere 
tmt  convenient  sines  they  are  too  matbemecical  and  too  abstract  for 
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formulations  and  aoaa  quanclclas  ara  not  directly  maaaurabla  in  the 
aquations.  In  ordar  to  obtain  useful  solutions » the  complicated  Integrals 
In  the  conservation  equations  have  to  be  simplified  by  the  proper  definition 
of  averaged  quantities.  The  final  form  of  this  simplified  equation  should 
agree  with  the  aquations  derived  from  the  tvo-phaae  fluid  dynamics  approach. 

The  gas-dynamic  approach  is  then  formed.  Such  an  approach  will  pro- 
vide a clear  definition  for  each  term  in  the  governing  equations  in  both 
phases.  Some  bulk  properties  of  the  fluid  defined  in  the  formulation  and 
the  gross  effects  of  heat  transfer  and  flow  resistance  can  be  measured 
experimentally . 


CHAPTER  II 


LITERATURE  SURVEY 

la  chi*  chapter  vt  provide  a literature  survey  on  three  different 
topics,  namely  (i)  a brief  survey  on  .stint  two-phase,  vsactlve  flow 
combustion  nod  els ; (il)  empirical  correlations  such  as  correlations  for 
heat  transfer,  flew  resistance  and  lntarfranular  stress;  (lli)  numerical 
methods  to  solve  this  aatheaatical  nodal. 

2.1  Survey  on  Related  Works 

The  objective  is  to  sismrite  the  past  axperies  as  of  various  in- 
vestigators working  in  this  fisld,  so  that  all  inport ant  physical  phenomena 
and  influencing  factors  era  included  in  the  theoretical  nodal. 

Before  going  into  several  recant  investigations  of  tr tux lent  burning 
of  granular  propellants,  so ns  published  work (37,  srj on  the  steedy-etste 
burning  of  slnilar  propellents  are  reviewed.  Steady-state  and  transient 
burning  studies  are  alnci  inseparable,  therefore,  eons  discussions  of 
current  steady-state  results  prove  to  be  worthwhile.  Kuo  at  al  (l,  39)  show 
that  gas -penetrative  burning  of  porous  propellents  or  explosives  under 
strong  confinement  is  Inherently  self-accelerating.  Kuo  and  Sunwrfield 
$7,  IS)  have  shown  that  under  suitable  physical  conditions,  the  internal 
pressure  distribution  generated  by  propellant  gasification  esn  produce  a 
constant  tats  of  gas  penetration  through  the  taburned  region,  and  nay 
eventually  give  a constant  spaed  of  combustion -wsvs  propagation. 

Several  investigators  (AO , Al,  42,  43)  have  Observed  experimentally,  soma 
evidence  of  steady  gas -penetrative  burning  in  porous  proep Hants.  Recent 
irtnsjits  mere  made  by  various  investigators  to  gain  a deeper  understanding 
of  the  physical  pftooomois  Is  1-nlting  add  burning  of  granular  ptopellsccs . 


In  eh*  consideration  of  ch*  *ff*ct  of  propellent  typ^  on  eh*  praseure- 
time-distance  variations,  Wer lick  (44)  conducted  several  experiment*  on 
various  artillery  systems  using  XACO  propellents  (the  cool-burning  pro- 
pellent which  reduces  barrel  erosion) , end  comparisons  were  made  with  the 
mo  propellents.  The  pressure- time  srtces  obtained  at  the  case  base 
showed  no  significant  difference  with  propellant  type.  The  pressure- 
time  history  at  the  projectile  base  le  later  found  to  be  radically  different 
from  those  recorded  at  the  cartridge  case  base.  The  Intensity  of  the 
traveling  pressure  front  In  the  cotabuation  chamber  was  found  to  be  strongly 
dependant  on  primer  venting,  porosity  distribution,  and  the  physical- 
chemical  characteristics  of  the  propellent  In  the  bed.  War] ink  firmly 
suggests  that  intarlo?  ballistic  systsas  studies  must  Include  the  con- 
sideration of  the  above  affects  to  minimize  the  pressure  irregularities. 

Soper  (lQ) performed  experiments  on  ignition  waves  in  gun  chambers 
loaded  with  NACO  propelle  r's.  Hit;  experimental  results  revealed  severe 
wave  action  during  Ignition:  a distinct  pressure  wave  traveled  down  and 

continued  oscillating  Inside  the  combustion  chamber.  A strong  prersur* 
spike  at  the  base  of  the  projectile  wee  recorded,  (the  typical  Intensity 
Is  2460.0  kg/ cm2.  Soper's  flesh  X-rey  apparatus  detected  e pronounced 
acceleration  of  propellent  particles  ahead  of  the  gas  pressure  rave.  This 
suggests  consideration  of  both  tha  motion  of  propellents  end  also  the  Intsr- 
greovler  forces  transmitted  la  compaction  beds  ere  equally  important. 

East  and  McClura(4)  set  up  a transient  interior  ballistic  model  by 
considering  the  two  phase  phenomena.  Their  mathematical  model  Is  not  quite 
complete  for  the  following  reasons.  They  assumed  that  tbs  initial  state  of 
tha  system  in  already  fully  Ignited.  This  assumption  defeats  the  major 
purpose*  of.  analysing  tha  offset  of  Ignition  waves  and  optimising  the  primer 
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design.  In  chair  |ai  noman tun  equation,  the  mmocu  of  riuct*^  ui a «u 
incor. act.  Tha  solid  phase  normal  atraaa  transmission  vaa  excluded  in  chair 
solid  phaaa  acmantua  aquation.  In  thair  gas  energy  aquation,  tbmy  neglected 
tha  hast  loaa  to  tha  propellant  from  tha  gaa  phasa.  Host  important,  tha 
•olid  phaaa  energy  aquation  was  excluded  free  their  formulation,  and  tha 
surfaca  temperature  of  tha  propallanta  vaa  not  treated. 

Gough  and  Smarts  (3)  mada  a good  attampt  in  formulating  a aathamatical 
modal  for  ignition  of  porous  propallanta  by  considering  tha  motion  of 
solid  particles  and  atraaa  transmission  in  tha  compacted  bad.  Thair  modal 
vaa  far  mors  sophisticated  than  the  East  and  McClure  modal.  However,  a 
few  errors  of  tha  medal  <ira  manlfastad  in  tha  fallowing.  Ergun's  correlation 
(43)  was  used  to  calculate  tha  totem  between  gaseous  and  solid  particles 
In  both  fluidised  cad  noa-fluidlxed  sections  of  the  granular  bad.  But 
*;^vs's  rralitlsa  is  ssly  valid  for  nos** fluidised  (fixed)  bed.  In  addition 
tc  the  iutasaal  thermal  energy,  lie  total  energy  addition  to  tha  gaa  phaaa 
doe  to  burning  should  include  bc*A  tha  chamlerl  energy,  and  tha  flow  work 
«aich  drives  tha  reacted  gaseous  mass  into  tha  space  occupied  by  other 
|;«u  in  .ee  control  volume.  Tha  flow  work  was  not  Included  in  tha  gas- 
phAt  a vwrgy  aqration  of  thair  modal.  Thair  boundary  condition  on  the 
projectile  and  after  tha  shot  start  mas  r^t  wall  specified.  Gough  and 
Ewarts'  numerical  solutions  indicated  samara  oscillations  an  the  boundaries. 
Thao#  spurious  waves,  game rated  by  numerical  Instabilities,  propagated  into 
the  interior  and  iaaccuraalaa  resulted.  Gough  (12,  id)  had  Improved  his 
modal  by  rv^exmainias  tha  csaaarvatlon  o^aatioaa  and  this  la  diseased  later 
in  this  seat ion. 

fisher  at  al  (2,  47,  4ft,  4ft)  stall  ad  propallamt  Igaltios  amd  combustion 
proa— aa  in  aswannl  different  artillery  07s tana;  tha  gaanrtxy  involved  vaa 
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quit*  different  from  « ona-dlaeaslonal  modal.  Th«y  considered  th*  combustion 
proc*3*o*  in  th*  ifniclon  cartridge , th*  flow  of  hot  r**ct*d  gases  and  pro- 
pellant radirlly  through  eh*  vent  hol«*  on  th*  vail  of  th*  ignition  cartridge, 
and  th*  combustion  of  th*  bag  propellant.  Although  th*  action  of  propellants 
la  consid*r*d  by  th*  conservation  aquations  for  th*  propellents,  the  modal* 
nr*  incomplete  and  their  predictions  have  not  he*m  verified  thoroughly  by 
experimental  data. 

frier  at  al  (5,  50,  99)  approached  the  tvo-phase  granular  propellent, 
burning  problem  in  a different  aennar.  They  took  the  continuous  mechanics 
approach  for  both  gaseous  and  solid  phase*.  They  treated  the  propellent 
particles  as  * fluid.  Eowever,  the  fluid  and  aolld  media  differ  distinctly 
b7  the  characteristic  differences  in  (tress  tensor.  In  e fluid  it  Is  pro- 
portional to  the  cate  of  angular  deformation,  vherea*  th*  stress  tensor 
in  a solid  medium  la  proportional  to  the  magnitude  of  deformation.  To 
analyse  the  ceusa  of  generation  of  the  pressure  spike  at  th*  base  of  th* 
prcjectil*  and  the  wave  action  In  tv*e  unhurried  propellent,  the  stress 
transmission  in  the  solid  material  muse  be  treated  differently  from  th* 
stress- to-strala-rnt*  relationship  for  fluids.  In  addition,  the  energy 
equation  for  the  solid  phase  was  written  in  terms  of  global  stored  energy. 

This  lumped  parameter  is  not  adequate  for  studying  the  flams  propagation 
sad  for  calculating  the  heat  transfer  from  the  hot  gs*  to  the  unbruaei  pro- 
pellents. The  conduction  process  inside  the  solid  propellent  la  a relatively 
•law  process  in  comparison  with  the  overall  trees  lent  characteristic  of 
fleam  spreading.  The  temperature  profile  within  the  solid  particles  must 
have  s sharp  gradient;  it  la  vagu*  to  use  e ltatped  paresmtsr  in  rha  solid 
phase  energy  equation.  The  investigators  realised  the  importance  of 
meltable  empiric*!  antclaelna  for  the  fluidised  esd  aon-fluldlsed  bed 
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and  they  collected  * brief  survey  on  existing  correlations  for  heat  transfer 
and  flow  resistance.  The  boundary  treatment  i*  inadequate  since  the 
reflection  condition*  are  u**d  at  th*  boundaries. 

tough  (12,  46)  in  hi*  sort  secant  study,  ha*  invsstigetsd  th*  interior 
ballistic*  of  guns  by  considering  the  microscopic  aquations  of  motion 
tiuwd  on  'uv.Mj  anti  explicit  av»*  epp'itca.  Si^  a^pr©*;!.  i*  piredietai 

upon  a combination  of  th*  nathodt  of  laderson  end  Jaekaon  (2?),  Whitaker  £}0)  *nd 
P«aton(l3l  A weighting  function  is  adopted  for  tha  definition  of  all  the 
average  flow  properties.  Overall,  this  proved  to  be  a useful  scheme  to 
dsslve  a set  of  conservation  equations  directly  from  the  aicroeccpic  point 
of  view  to  give  the  flow  properties  in  the  microscopic  sense,  hut,  this 
approach  is  too  aathcaatical  and  too  abstract  to  give  tangible  physical 
meanings  to  each  term  in  tha  conservation  aquations.  Cough's,  modal  still 
uses  the  Erps?  equation  f a?  fluidised  sad  packed  bads.  It  has  not  bean 
verified  by  some  of  the  rigorous,  one-dimensional  experiments . With  s 
continuous  impi  ovement  in  th*  boundary  point  treatment  and  accurat*  forcing 
function  specification,  it  is  expected  that  Gough's  model  should  give  dose 
predictions  to  th*  results  of  our  model. 
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2.2  Survar  of  EgglricjU.  Correlations  mad  Their  Limitations 


Ir  coo 'true ting  a complete  thaoratical  nodal,  some  information 
baaidaa  tha  mojor  ernservation  aquations  ara  neceesary.  These  anpirical 
correlations  aarve  to  dacamlna  haat  transfer  in  tha  granular  bai,  drag 
forces  between  partldaa  and  gaaaa,  and  atraaa  trananittad  threugh  closely 
packad  granular  prop a 11 ant a . 

An  axtanalva  literature  survey  was  conduct ad  to  collact  availabla 
corralationa  on  haat  tranafar  and  drag  data  for  non-fluldixed  and  fluidized 
bada.  Tha  chaalcal  angluaarlng  industry  ganaratad  tha  majority  of  tha 
availabla  corralationa . A largo  quantity  of  fluidization  studlaa  hava  baen 
contalnad  In  hundrada  of  raporta  and  papars,  and  alao  in  taxtboc&a  by 
Zanx  and  Othaar  (51)  £ava  (52)  £enu  and  Lavanaplal  (53)  ,Soo  (54)  ,and  Davidson 
and  Sarrlaon  (5$  . Unfortunately,  empirical  corralationa  that  covar  a wida 
ranga  of  Raynolda  numbers  and  poroaitlaa,  partida  sixa  and  alao  lnduda 
tha  affact  of  partida  (taaif lcatlon , ara  unavallabla.  Krlar  at  al(s) 
ancouatarad  tha  a ana  pro bias  in  thair  racant  raport.  Aa  a raault,  tha 
empirical  corralationa  for  flow  ovar  inart  spharaa  in  non-fluidizad  and 
fluidized  bads  ara  employed  in  various  combustion  modal a,  and  even  extra- 
polation to  larger  ranga  of  Reynolds  numbers  and  poroaitlaa  to  continue 
with  the  calculations  on  a computer. 

Since  no  suitable  drag  correlation  is  available  for  a vide  ranga  of 
Reynolds  number  and  porosity,  wa  performed  a large  number  of  cold  flow 
experiment*  and  obtainad  a useful  correlation  between  tha  pressure  drop , Raynolda 
number  and  porosity  for  non-fluidlzad  bada  of  VC  870  propellants  (7$. 

Tha  Raynolda  number  has  bean  extended  to  be  six  timss  higher  than  that  in 
Ergua'a  aquation.  For  fluidised  regions,  Andsrsaon's  Correlation  w6)  haa 
bean  found  to  bm  a suitable  anpirical  formula,  sines  it  is  valid  for 
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porosity  from  0.45  to  1.0.  However,  it  must  b«  noted  that  the  Feynolds 
number  verlee  only  froti  0.003  to  2000  in  hie  correlation.  It  must  be 
used  with  caution. 

For  convective  heat  transfer  calculations, Denton’ s (57)  formula  is  found 
to  be  adequate  for  the  n on-fluidized  region.  It  is  valid  for  Re^  from 
500  to  50,000  and  porosity  around  0.37.  Rowe-Claxton ’ s correlation  (5?)  is 
found  to  be  the  best  one  for  describing  the  convective  heat  transfer  in 
fluidised  regions.  This  expression  is  valid  for  Re^  froa  10  to  IQ7,  and 
oorosity  froa  0.26  to  0.632. 

A constitutive  lav  is  needed  to  specify  the  normal  stress  transmission 
through  packed  granular  particle* . When  particles  are  dispersed,  the  inter- 
granular stress,  T is  aero.  When  the  particles  are  packed,  T may  be  large 
P P 

due  to  the  physical  contact  between  the  particles  under  steep  pressure 
gradients . 

After  e literature  survey  on  intergranular  stress  correlations , it  is 
found  that  not  much  information  is  available.  Soper  0-Q)  observed  that  the 
force  and  displacement  of  a piston  acting  on  a granular  bed  of  NACO  pro- 
pellents, were  llnaarily  relatad  vxthln  the  range  of  0 to  351  compression. 
The  speed  of  sound  in  granular  propellent  was  dsduced  from  hie  experimental 
data.  Besides  Soper's  vork,  the  dependence  of  particulate  stress  tensor 
upon  porosity  of  the  mixture  was  given  by  Ho.et^S)  and  Gough 0-9*  Further 

discussion  on  T Is  gives  in  Appendix  6. 

P 
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2.3  Survy  on  Numerical  Studies 


After  the  theoretical  formulation  is  completed  and  the  nature  of  the 
governing  equatlona  is  studied,  a literature  survey  in  selecting  a stable 
and  fast  convergent  numerical  scheme  has  been  conducted. 

Obviously,  the  theoretical  model  is  a system  of  coupled,  non-linear, 
inhomogeneous  partial  differential  equations.  Since  the  eigenvalues  are 
all  real,  the  system  is  hyperbolic.  In  order  to  select  a suitable 
numerical  scheme,  two  important  criteria  must  be  noted  in  making  such  a 
decision.  They  are  the  stability  and  convergence  of  the  scheme.  Papers, 
reports  as  well  as  textbooks,  have  numerous  types  of  finite-difference 
schemes  to  handle  hyperbolic  partial  differential  equations.  Moretti(6Q) 
make a a good  comparison  of  the  differences  of  some  r merical  schemes.  The 
stability  criteria  has  been  studied  for  each  scheme.  In  practice,  a 
rigorous  stability  analysis  of  the  sat  of  coupled  non-linear  partial 
differential  equations  cannot  be  performed  due  to  the  complexity  and  non- 
linearity of  the  equations.  For  the  study  of  non-linear  systems  of  equations, 
Hicks (6^  suggests  disregarding  the  study  of  stability  criteria  and  proceeds 
directly  to  the  heart  of  the  problem,  namely,  the  convergence.  This  was 
also  discussed  by  Lax  and  Richtmyer  (62)*  Basically,  the  FDE  (finite-difference 
aquation)  solution  must  approach  the  PDE  (partial  differential  equation) 
solution,  and  stability  criteria  then  becomes  of  seconder,  interest.  As 
a result,  several  convergence  tests  were  suggested  to  test  the  numerical 
solution’s  convergence  to  the  unique  solutions  of  the  system  of  partial 
differential  equations. 

None  of  the  criteria  or  analyses  srs  adequate  for  practical  computations. 
In  actual  fluid  flow  problems,  the  stability  restrictions  are  applied 
locally.  The  typical  and  prudent  practice  is  to  use  some  percentage,  usually 
80%  or  90%  of  as  Mialytieally- indicated  maximum  time  step.  When  transients 
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are  large,  smeller  percentages  are  needed  to  stabilise  tha  solution.  For 
example , Terrence  (63)hes  dons  sobs  computations  on  natural  convection  by 
using  a s&allsr  At  to  achieve  computational  stability. 

The  shortcomings  of  the  above  approach  ere  not  very  apparent.  In- 
vestigators such  as  Phillipa(&4)»  Rishtayer  (65).  Hirt  (66)  and  Gourlay  and  Morris 
(67)  have  reported  instabilities  from  nonlinearity  or  from  non-constant 
coefficients.  Lilly(68)has  reported  the  phenomena  of  tine  splinting  of 
solutions  'rfhlch  is  not  precisely  an  instability  in  the  sense  of  producing 
unbounded  solutions,  actually  an  instability  due  to  interactive  convergence. 
Therefore,  it  is  important  to  realise  that  no  distinct  difference  exists 
between  a "tine"  instability  and  a poor  convergence. 

Overall,  four  methods  are  studied  in  detail  to  evaluate  the  stability 
of  tbs  numerical  s chare,  namely,  the  Courant-Friedrlcka-Lswy  (C.F.L. ) (69) 
conditions,  tha  Inargy  Method  (7 Q).  the  von  Neuiaxm  method  and  tha  Heuristic 
Stability  Theory  of  Hirt  (6$ 

In  the  selection  of  numerical  schema  for  tha  hyperbolic  governing 
equations,  several  different  schemas  are  considered:  the  method  of  Courant, 

Isaacson  and  Rees  (7$  the  Lelevler  scheme  (7C^  the  two  step  Lax-Wendrof £ 
scheme  (65;  ,the  * «p  Frog  scheme  (2 ail  the  two-step  Richtmyer  scheme  (oq),  the 
generalised  implicit  scheme  (7$,  the  method  of  characteristics  (74)  and 
tha  method  of  linos  (78).  From  the  previous  investigators'  experience  we 
finally  decided  upon  a modified  two-step  Richtmyer  schema.  This  will  be 
discussed  in  detail  in  Chapter  17. 
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CHAPTER  III 


THEORETICAL  FORMULATION  AND  ANALYSIS  OF  TRANSIENT 
COMBUSTION  PROCESSES  IN  GRANULAR  BEDS 

We  begin  this  chapter  with  tha  description  of  the  physical 
processes  Shat  occur  in  an  actual  firing.  Then,  after  a brief 
discussion  of  sons  basic  assumptions  in  this  nodal,  we  formulated 
our  taathenatical  nodal  by  deadrlbing  tha  physical  meaning  of  each 
tern  in  the  conservation  aquations  of  mass,  momentum  and  energy  of 
the  gas  phasa;  and  tha  conservation  aquations  of  mesa  and  momentum 
of  the  particle  phasa.  Tha  energy  equation  of  tha  particla  phase  is 
treated  separately.  Tha  conaarvation  aquations  derivsd  in  section 
3.3  ere  vary  general  in  nature.  In  order  to  arrive  at  e set  of 
equations  that  car  ha  solved  economically  on  tha  computer,  various 
amplifications  are  made  in  section  3.4. 

A set  of  governing  equations  it  tha  entrance  region  arn  derived 
to  facilitate  the  calculations  at  the  laft  boundary.  Tha  caturs  oZ 
the  governing  equations  are  determined  in  section  3.6  end  detailed 
mathematical  steps  are  shown  in  Appendix  3 and  4.  Initial  and  boundary 
conditions  ara  apeciflad  in  section  3.7,  further  discussions  on  the 
numerical  treataant  of  tha  boundary  conditions  are  given  later. 

In  order  to  close  this  problem,  associated  empirical  correlations 
and  limitations , such  ms 

A.  Conatitutiva  law  for  tha  total  stress 

B.  Drag  correlation 

C.  Heat  transfer  correlation 

D.  Burning  lata  law 

and  the  ignition  criterion  for  the  solid  propellent  are  separately 
discussed  la  thv  remaining  of  tha  chspr.tr. 
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3.1  Description  of  Fhvslccl  Prop 


Its  preparation  for  an  experimental  firing,  spherical  granular  (VC  870) 
propellent*  art  loaded  into  a thick-welled  cylindrical  steel  chamber  with 
15.24  cm  is  length  end  0.777  cm  I.D.  (Fig.  1).  At  the  left  end  of  the 
cylindrical  chamber  there  is  a gaseous  pyrogen  Igniter  which  Introduces 
hot  gas  into  the  gr&nular  bed  by  forcing  it  through  a multiple  perforated 
nozzle  to  achieve  a uniform  stream.  At  the  right  end  of  the  cylindrical 
chamber , there  is  e 0.81  am  stainless-steel  shear  disc  ’’burst  diaphragm", 
(this  provides  a fixed  boundary).  Ae  the  igniter  gee  is  turned  on,  the  hot 
igniter  geuas  generated  cause  e compaction  of  the  granular  bed  near  the 
entrance  region.  This  also  heats  up  the  nearby  granular  propellents 
to  ignition  condition.  The  burned  propellents  give  off  more  hot  genes, 
tbasa  gaseous  products  are  then  pushed  forward  by  the  pressure  gradient 
to  ignite  sore  propellents.  Thus,  a sce«p  pressure  gradient  i*  created 
Inside  the  combustion  chamber.  The  rate  of  presevrisation  increases 
further  in  the  downstream  direction,  resulting  in  a rapid  acceleration  of 
the  flame  front . The  shear  disc  begins  to  rupture  when  the  preei  ure  in 
front  of  it  reaches  e critical  value.  After  the  sheer  disc  bursts,  the 
combustion  chamber  pressure  decreases  abruptly  and  terminates  the  event. 


i. 2 Baltic  Assumptions 


There  are  3 basic  assumptions  pertinent  to  this  present  study: 

1.  The  flow  in  the  granular  bed  is  considered  to  be  quasi-one- 
diaensional.  The  effect  of  fluctuations  in  the  f lane  front 
for  tvo-disansional  or  three-dimensional  model*  may  be 
obtained  vLta  further  aodifi cations.  In  this  nodal  the 
changes  of  flow  properties  in  the  transverse  direction  are 
considered  in  the  esplrical  correlations  for  heat  transfer 
are  flow  resistance*. 

2.  All  chemical  reactions  occur  on  the  propellent  surface, 

in  the  conbustion  zone,  a thin  region  that  is  very  small  in 
comparison  with  initial  particle  diameter. 

3.  The  dynamic  burning  effect  ia  negligible  i.e.  the  rate  processes 
at  tha  pro ep lien t surface  are  quasl-stealy  in  the  sense  that 
the  characteristic  times  associated  with  the  gaseous  llama  and 

A 

preheated  propellent  an  short  compared  to  that  of  the  pressure 


transient  variations 


3.3  Theoretical  Modal 


In  order  to  study  the  detailed  structure  of  the  transient  gas 
dynamic  behavior  of  hot-gas  penetration,  flame  propagation,  chamber 
pressurization  and  combustion  processes  in  the  granular  propellent 
bed,  a set  of  governing  equations  describing  the  changes  of  mass, 
momentum,  and  energy  for  the  gas  phase  and  the  changes  of  mass, 
momentum  for  the  solid  phase  are  derived.  This  set  of  governing 
equations  is  expressed  in  a quasi -one-dimensional  form.  The  equations 
are  approeched  by  considering  the  belance  of  fluxes  over  a control 
volume  small  enough  to  give  the  desired  spatial  distributions  in  the 
complete  system,  yet  large  enough  to  contain  many  solid  particles,  so 
that  the  averaged  particle  velocity  and  fraction  porosity  are  meaning- 
ful. 

The  control  volume  for  the  gas-phase  is  the  portion  of  void  volume 
occupied  by  the  gas  phase  in  a small  eleuentary  volume  (A  Ax  )»  while 
the  remaining  portion  occupied  by  the  particles  la  considered  to  be  the 
control  volume  for  the  particle  phase  (Fig.  2). 

A Conservation  of  Hass  Equations 

In  deriving  the  gas  phase  conservation  equation  the 
fractional  porosity  is  defined  as 


a-  to1u**  - 1 adfcoS 

total  volume  *1  A AX 

- | - 


(3.3.1) 


Where  n is  the  number  of  pellets  in  the  cylindrical  volume 
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A A x,  and  r Is  tha  pellet  radius  in  cha  control  volume.  n is  the 
P 

number  density  of  spherical  pallets  in  the  granular  propellant  bed. 
The  specific,  surface  area  which  represents  the  total  surface  area  per 
unit  spatial  volume  Is  given  by 


ulmuLi 

A ax 


rP 


(3.3.2) 


The  gas-phaae  mess  equation  Is 


M£&2  + .£IP“J01 As  p-r>  (3.3.3) 

d t, 


The  rate  cf  gaseous 
aass  accumulation 
In  the  control 
volume  occupied 
by  gases 


The  rate  of 
Increase  of 
gaseous  mass 

by  convection 


The  rate  of 

gaseous  mass 

addition  due  co 
gasification  of 
solid  particles 


Tha  parti  da-phase  mass  equation  Is 


aCClrdM+  m&eezd- 

at  ax 


-A&n, 


(3.3.4) 


The  rate  off 
particle  aass 
accumulation  In 
tha  control  volume 
occupied  by  par- 
ticles 


Tha  rate  of  Tha  rate  of  partlde 

increase  of  aass  reduction  due 

particle  aass  to  gasification  of 

by  convection  solid  partidas 


In  this  analysis,  the  gas  density  p la  based  on  the  gaseous  mass 

divided  by  tha  void  volume  occupied  by  tha  gsaes.  The  propellent 

density,  o is  assumed  to  be  constant. 

P 


B Tha  Conf •rvation  cf  Moaactua  Eguatlcna 


The  gas-phase  »a**a tun  aquation  la 


a «t>P“s) 
at 


3 (OP**') 

* a V 


+ 3 


MJ&E2- 

ax 


Tbs  rata  of 
increase  of  gaseous 
■oMotuB  in  tha 
control  voluaa 


Tha  rata  cf 
■oatntua  flux 
incraasa  in 
tha  control 
voluaa  by 
eonvaction 


Pros aura  gradient 
fores  acting  in 
tha  x-dlrectlon 


a 9(T„g>)  _ 

3 ax 


AsrbPp*f  - 9A,Dt 


Viscous  fores 
dua  to  normal 

stresses 


Tha  rata  of  In-  The  total  drag  fores 
crease  of  linear  between  the  gas  and 
aoaentun  of  gas  particle  phases 
phase  due  to 
the  notion  of 
gasifying  par- 
ticles in  the 
aoaentua  direction 


~ g (3.3.5) 

Friction  force 
acting  oa  the  gae 
by  tha  charter 
well 


» 


Th*  part.lcla-phaaa  momantu a • qua  t Ion  la 


aCti-0)P.Up]  + afo-0?P,V] 
at  ax 


" ax 


Tha  rata  of 
lncreaaa  of 
partlcla  mo- 
mantus  la  the 
control  values 


The  rata  of 
momentum  flux 
lncraaaa  is 
tha  control 
volume  by 
convection 


Intergranular  force 
acting  on  particle 
cross-sectional 
arena  on  the  boundary 
of  the  control  volume 


— Asr;p,up  + gAsDt 


A (3.3.6) 


The  rate  of 

decreaae  of 

particle  mo- 
mentum due  to 
gaelficatlon  of 
particles  into 

gas  phaae 


Tha  total  drag  Friction  force  acting 
force  between  on  the  particle  by 
the  gas  and  the  chamber  wall 

particle  phase 


It  is  important  to  not*  that  is  the  total  drag  force  between 
ft i gas  and  particle  phases.  It  is  equal  to  the  sum  of  the  drag  due 
-am  presence  of  relative  velocity  between  the  gas  and  particle 
prases  and  the  drag  due  to  the  porosity  gradient,  l.e. 

Dt-  Dv+  Dp  - Dv“  ^ C2.3. 7) 

The  expression  of  will  be  discussed  in  a later  part  of  tha  theoratlcal 
mode)  . The  advantage  of  writing  the  above  two  momentum  equations  In 
texv/s  ef  D,  la  to  facilltats  tha  solution  of  tha  hyperbolic  problem 
which  la  not  always  totally  hyparbollc.  Tha  alganvaluas  may  net  ba  all 
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distinct  whan  tha  granular  bad  bacon**  disparate!.  Tu*  aecond  ad- 
vantaga  la  to  facilitat*  tha  transformation  of  tha  governing  aquation* 
into  tha  charactariatic  form  and  to  daarly  indicat*  tha  conpatibillty 
relations  for  tha  flow  propartlaa  at  tha  boundarlaa  of  tha  granular 
bud. 


C Tha  Conservation  of  Energy  Equation 
Tha  gaa-phaaa  anargy  aquation  la 


2L&EM  + ZSS££Va£L  + l afPuggj 

at  ax  j ax 


Tha  rata  of 
lncraaaa  of  total 
storad  anargy  In 
control  voluna 


Tha  rata  of  anargy  Tha  rata  of  vork 
flux  Incraaaa  by  dona  on  tha  gaa 

convection  phaaa  by  prasaura 

fore* 


_ a<  qgn 
ax 


-A.htlT-T") 


Tha  rata  of  baat 
conduction  through 
tha  control  aurfaco* 


Work  dona  par 
unit  tl mi  by 
tha  total  drag 
£orca  acting 
on  partial** 


Baat  losa  to 
partlclaa  par 
unit  tin* 


+ Pp rb  (! ^charn  ♦ ) 

Sh»  gaga  cftatal  anargy 
input  hra  to  tha  cowfcuatlen 
•f  partial** 


, i XT«u«t» 

T .ax 

Tha  rata  of  mirk  dona 
hy  via con*  normal  atraaa 
in  tha  gaa  phaaa 
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(3.3.8) 


P 3(0) 

J at 

Tha  rata  of  praasura  work 
for  tha  dilatation  of  tha 
gaaaoua  control  volume 


-Q 

wvw 

Tha  rata  of  heat  loaa 
to  tha  combustor  wall 
par  unit  voluaa 


Tha  empirical  corralation  for  h£  will  ba  shown  latar . Tha 
paraaatar  J is  tha  haat-to-work  unit  conversion  factor.  Tha 
conservation  aquations  in  Eularian  font  ara  shown  in  Appendix  1. 

Tha  abova  sat  of  govarning  aquations  are  darlvad  in  tha  fora  of 
coup led « non- linear,  inhomogeneous  partial  diffarantial  aquations. 
Further  simplifications  ara  required  to  solve  than  economically. 

D Solid-Phase  Heat  Conduction  Equation 

▲ particle-phase  energy  equation  similar  to  Eq.  (3.3.8)  can 
ba  darlvad  in  tors*  of  tha  average  bulk  temperature  of  a typical 
solid  particle  in  tha  particle-phase  control  volume  considered. 
However,  since  the  transient  interval  of  the  physical  procaaa  is  so 
short,  tha  ;aoperature  pro  fils  ins  id  a tha  particles  has  a vary  ataap 
gradiant.  Although  the  average  bulk  temperature  of  a solid  particle 
can  be  obtained,  it  would  nut  be  useful  for  determining  the  Ignition 
condition,  the  speed  of  flame  front  propagation,  or  the  rate  of  heat 
transfer  from  tha  gaa  to  solid  phase.  Zastaad  of  traatlng  tha  solid 
phase  anargy  equation  like  the  mass  and  a omen  turn  equations  for  tha 
particles,  the  transient  heat  conduction  equation  la  considered,  so 
that  tha  propellant  surface  temperature  and  the  temperature  gradient 
cam  be  accurately  obtained. 
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Following  th*  motion  of  a given  pa.:  tide,  the  heat  equation 
written  in  the  Lagrangian  time  derivative  and  spherical  coordinate 
of  the  particle  is  (see  Fig.  25) 


tm a)m  sL:  G'lTTr) 
iDt % r ir* 


(3.3.9) 


The  initial  condition  for  equation  (3.3.9)  is: 


t-0  : TP(0,  D-To 


(3.3.10) 


The  boundary  conditions  before  the  thermal  wave  has  penetrated 
to  the  center  of  the  sphere  are: 

r— o : ^£(t,Q)  - 0 (3.3.1D 


r-  r* : JfftL[rm-Tf3(t)J  (3 

where  the  total  heat  transfer  coefficient,  fc,,,  is  the  in  of  the  con- 
vection and  radiation  coefficients;  l.a. 

(3.3.13) 

mere  the  gas  eniaaivlty  is  taksn  to  be  unity. 

The  heat-up  equation  and  its  initial  and  boundary  conditions  are 
rec ea tad  by  the  integral  method  as  shown  in  Appendix  2 to  yield  a 
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first  order  ordinary  differential  equation  describing  the  Increase 
of  propellent  surface  temperature  with  respect  to  tine,  l.e. 


where  6 is  the  thermal  wave  penetration  depth  in  a spherical  particle, 
measured  fron  the  particle  surface;  £ as  a function  of  time  is  given 
by: 


g/x, arp.rw-tKT.J 

J-(Tbs(W-T.J+  dkj^[T(t>-Tre(t)J  «u.u> 


Similarly,  an  ordinary  differencial  equation  describing  the  changes 
of  after  the  thermal  wave  has  penetrated  the  full  radius  of  a 
spherical  particle  is  obtained.  Since  the  transient  interval  of  the 
physical  process  is  only  a few  nilllaecords , the  particle  surface 
reaches  the  ignition  condition  long  before  the  thermal  wave  penetrates 
to  the  center  of  the  particle;  the  equation  describing  the  changes  of 
after  tha  thermal  wave  penetration  is  not  given  hers  but  will  be 
shown  in  Appendix  2.  After  Tpg  is  solved,  a simplified  ignition 
temperature  criterion  is  uae^  to  determine  the  burning  condition  of 
particles  along  the  propellent  bed.  The  ignition  criterion  ie  given  in 
section  3.9. 
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3. A Simplification* 


3*for*  th«  necessary  boundary  conditiuna  art  considered,  tha 
governing  aquation*  for  tha  gaa  and  partlda  phases  are  almpllfiad 
furtbei  after  an  order  of  magnitude  analysis.  The  highex-ordar 
tarn*  neglected  ere:  (a)  the  viscous  normal  stress  in  the  gaa-phaat 

momentum  equation,  (b)  the  shear  force  at  the  combustor  wall  for 
both  gases  end  particles  in  their  momentum  equations,  (c)  the  gas- 
phase  heat  conduction  tern,  (d)  the  work  done  by  the  viarou*  normal 
stress  in  the  ga s-phaae  energy  equation,  (e)  the  boat  loss  0 to  the 
chamber  wall  in  the  extremely  short  transient  combustion  experiments, 
(f)  the  rate  of  pressure  work  for  the  dilatation  of  tha  gaseous 
control  volume  in  the  gas-phase  energy  equation,  end  (g)  the  rete  of 
change  of  the  total  heat  transfer  coefficient  in  the  calculation  of 
propellent  surface  temperature,  (h)  body  force,  (1)  dynamic  burning 
effect  cm  the  burning  rete  calculations  (j)  temperature  dapendence 
of  constant  pressure  specific  heat,  (k)  heat  dissipation,  (1)  con- 
stant cross-section'll  eras,. 

After  all  these  assumption* , a simplified  s*t  of  governing 
equations  will  he  rewritten  in  th«  following. 

The  ges-phuoe  aquation  hecomss , 


3m+2&m.Asp>n. 

dt  o 


(3.4.1) 


Tha  particle-phase  mass  equation  becomes , 


, 9U\-mu pj 

Zt  2X 


—AsfrCb 


(3.4.2) 
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Th«  gaa-phaae  momentum  equation  becomes. 


3(0PU»>,  + 4 2&PI 

3t  3X  3X 


-A*  n>fi>Up~gA*Dt  ( 3 . & . 3) 


Tha  particle-phase  momentum  equation  be cone c , 

afft-rtfWJ,  aftWjWJ  naia-0)T,] 
at  <3*  * ax 


--AiFeP^p  +gAsDt 


Tha  gas-phase  energy  aqua don  becomes , 


3(gP/^  + t i aiPtekL. 

at  + ax  f T ax 

, (3.4.5) 

— - Asht  (T-  T„)  + A %?p  ft  ^ehem^ 


The  particle-phase  heat-up  aquation  becomes. 


Tha  above  equations  (3.4.1)  - (3.4.6)  are  the  simplified  form 
r> { tha  governing  equations  (3.3.3)  - (3.3.6),  (3.3.8)  and  '3.3.14). 
These  equations  are  coded  directly  into  the  computer  program. 
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3 . 5 Governing  Eqoaricns  An  ch«  La ft  Control _VoJLuay. 

Tbs  tut  equation  of  the  gas  phatt  in  the  left  control  volume 
(Fig.  3)  is 


[ihijni't)-  rr,m(  t,  Xj] + Ase.  ^ Ph 


(3.5.1) 


wfcbre  is  the  ignitnr  uses  flow  rete  into  the  combust  ion  chamber 

end  4^  ( 

of  the  granular  bad,  it  is  represented  by 


t,  x^)  is  the  arss  flow  rete  of  gas  into  the  internal  ben  and ary 


min  (t.XJ-Plt.XjVj  (t.X.)A 0H.XJ 


(3.5.2) 


The  lease  equation  of  the  particle  phase  la 


<3.5.3^ 


Sene  as  above  ell  the  parsec teru  with  a subscript  c ere  evaluated 
ir  the  entrance  section. 

The  aonentvo  equation  **f.  the  gas  these  is 


+Adctr^^sw#-g^ 


U.5.4) 


where  u , the  Igniter  gee  velocity  is  given  'ey 

iiB  im^j 

(3.5.3) 
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Similarly  the  momentum  equation  of  particle  phaee  it  derived  in 
the  differential  form, 

dCCI-feyk r fl-C i 

dt  1-  E V 


+ASt(ga-rk£PreUfe) 


(3.5.6) 


For  the  control  volume  conaidared  in  Fig.  3r  the  energy  equation 
of  the  gaa  phaae  in  the  differential  form  it 


-pu9A0CCf7*^-)\J 

+ AiC  l-  -ht.CT.-W 

% 


(3.5.7) 


where  Efi  la  the  total  atored  energy , the 
aoergr  which  la  repreeented  by 

£t  -C» Tt+  wjj- 


aua  of  thermal  end  kinetic 


(3.5.2) 


h. Is  tha  enthalpy  of  th«  ignitar  gas  at  ths  tanparatura  of  the 

lgn 

hot  ignitar  gas  and  hch#m  is  tha  anthaipy  of  tha  propallant  gas  at 
flans  taaparatura. 

Tha  last  aquation  balng  consldarad  Is  tha  solid-phasa  hast 
aquation  In  two  forma 

(a)  Tina  Parlod  1 


r(6r*-Sc).  h,Al  m 

l r,.  + k,  j 


(b)  Tina  Parlod  2 


(3.5.10) 


In  tha  laft  control  vo liana,  thara  ara  a total  of  six  unknowns: 

4 , p*.  u , u , S , T 
« V gc'  Pc’  c’  psc 

Proa  Bq.  (3.5.3),  $ Is  axplldtly  daflnad. 

c 

p Is  Chan  daflnad  by  substituting  Eq.  (3.5.3)  Into  Eq.  (3.5.1) 
c 

It  bacosMs 
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(3.5.11) 


dPc 

dt 


I f ( WiinttyrhjnCt)) 

•*l — ax: — 


♦>ww.  #] 


Slailarly  fro®  £q.  (3.5.4)  w«  obtain, 

^ + Ajt(  Pkt  Pft  u*  -get) 

~Pc.UJt  '&***<.$?}  (3.5.12) 


By  th«  tarn  tokan,  Eq.  (3.5.7)  baeoaaa 


dEc_  1 fir nri. /+VK—  ) Du-  AOVCJJ*.  «ilil  1 

« JJRl>SD.  ' s*r  t ' ’ i*T'WJ 

+aJt  ^ -httn-ynj+wit w,^;J 


(3.5.13) 


Than  taaparatura  can  ba  a-valuatad  by  raarranglag  Eq.  (3.5.13), 

•lag 


Te  J^y”]  (3.5.14) 

tha  aquation  of  atata  fro®  Eq.  (3  1.1)  aftar  raarzaaglag  va  gat. 
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Pt_  Rlii 
(-PTb) 


(3.5.15) 


la  summery  <t  , u , u , T , P , T _ art  found  from  Eq.  (3.5.3), 

Pc  P c 

Eq.  (3.5.4),  Eq.  (3.5.6),  Eq.  (3.5.8),  Eq.  (3.5.9)  or  (3.5.10)  and 
Eq.  (3.5.15)  respectively.  All  these  parameters  are  us  ad  to  solve  for 
tha  lafe  boundary  condition*. 

Anochar  sat  of  govamlng  aquations  ara  darivad  at  the  laft  control 
volume  whan  there  ara  no  particles  existing  in  tha  entrance  region. 

Whan  there  ^re  no  particles  in  tha  region,  that  Implies  that 


0c- 1 


(3.5.16) 


<JLPc-0 


(3.5.17) 


Tha  maea  aquation  of  gas  phase  reduces  to 


dft  _ (miart-rn<n) 
« AX. 


(3.5.18) 


cad  tha  momentum  aquation  of!  tha  gas  phase  becomes 


4i 


As  s result , energy  equation  btcosu, 


+ -pa5A0(CrT-  jrfo)LL]/AXL 

-ec<L&-\  u,5,; 

c dt  J 


•nd  Tp#  calculation  is  slaplar,  by  ths  tine  all  the  particles  are 
pushed  out  of  the  entrance  region,  the  ignition  temperature  of  the 
particles  hsve  already  been  reached, 

Twc  *“  T»jn  (3.5.2 

The  above  equations,  (3.5.16)  - (3.5.21)  tfilx  bu  used  when  the 
left  control  vo  lusts  contains  no  particles. 


3.6  Slganvalua  Data ralattion  of  tha  Systsm 

To  study  tha  oatura  of  tha  govarnlng  aquations,  tbs  slgsnvsluss 
of  ths  systsa  wist  bs  datsrminad.  By  substituting  ths  governing 
aquations  In  ths  Eulsrias  fora  Into  s matrix  form 


• 

U3 

* 

-Us 

0 

<1 

o 

° 0P 

' 

Ufl 

' I 
I, 

T 

0 

0 0 

T 

It 
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0 

JL. 

ax 

P 

+ 

Ij 

O 

0 

0 

T&r 

Up 

u 

0 

o 

0 

0 

(£-0)  - Un 

0 

I* 

t 

i 4 

• 4 

(3.6.1) 


whsro  I „ I , I , I , I ara  ths  inhoaogansous  earms  of  tha  govaming 

1 2 J % 1 


aquations  and  ara  daflnad  as 

l,-fy[P'r*(Up-^)-3D>] 


♦ sijf^vflvw  > Pf  ofw  W * rfr^ 


T P> 


hmA*r> 


(3.6.2) 


Nov  w«  define  the  aetrix  M as: 


•“!  » f 

wz ; -**•  0 0 0 

M“  0 -US  0 (3.6.3) 

0 0 0 -Ur  7& 

0 0 0 o-tf>  -up 

. I - 

by  Mttin|  | M • X 1 | » 0 «i  can  obtain  tbs  elgearaluas . Detailed 
aatheaatieal  steps  era  shove  in  Appendix  4 . 

The  six  roots  found  are  aa  follows; 

A,  **—  tty +C$) 

At (Uf 'Cfl) 

A>  “-tfj 

A^--(^P^)  (3.6.4) 

J\$  — (U+~C) 

A.  ■*"  Up  (obtained  directly  by  following  the 
notion  of  a given  solid  particle) 

Since  all  these  six  roete  are  distinct  real  Gushers  for  non- fluidised 
oenditian,  the  eyetae  of  governing  ecus  titans  is  totally  hyperbolic  (97). 

Tbs  slope  of  the  characteristic  lines  for  the  governing  equations 
are  iftfsmrlired  by  tbs  local  gas  velocity  and  the  spaed  of  noad.  Its 
tmiMristk  dinc4«M  are  describee  by  9 


(If),""  " u* 


(3.6.5) 


A,  - Uj-C3 

drl—  A>  “ 
(KHf~A+  “u>'+c 


/dx  \ * 

a 5 


CdFT 


1 • /• 

up  — 


(3.6.6) 


(3.6.7) 


(3.6.8) 


(3.6.9) 


(3.6.10) 


In  the  nbcoele  cm,  Eq.  (3.6.5),  Eq.  (3,6.6)  and  Eq.  (3.6.7) 
npnant  the  right-running , la  ft- running  and  gaseous-path  Hach  linns, 
reapectlTeljr , In  fib*  gas  phase.  Iq.  (3.6.8),  Eq.  (3.6.9)  and  Eq.  (3.6.10) 
raprsaant  tha  right-running,  left-naming  and  particle-path  Mach  laoaa 
la  the  aalid  nhm 

naadawnl no  ell  Eq.  (3.6a)  lass  Its  characteristic  fern  Is 
fauad  by  calm*  the  ajgaa-dtra ctlena  ef  atari*  V as  a sat  of  land 
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variable  coordinates  in  tha  (t,  x)  planaa (J4  f 97),  cha  detailed 
mathematical  a tape  are  given  in  Appendix  3. 

Tha  characterietic  aquation  along  the  right-end  left-running 
Mach  linae  in  tha  gas  phase  are, 

(dPli*  Kiptd  ap)xfK  0 0)imKtm(eit)!t  (3.6.11) 
and 

4 K ?(dPh  +l<UdUp)*  (3.6.12) 


Tha  characteristic  aquation  along  tha  gaseous-path  Mach  line  in 

the  gas  phase  is 


. — 3 i/X/  1 «»>\  1/1/  ^1. . % , _!/•*  tAA\ 

f\T  VO  (3.6.13) 


Again,  the  eharactarlatic  equation  along  the  right-running , left- 
running Mach  lines  in  the  solid  phase  are, 

(<*«*)»  4'K*(^PJK  )&  (3.6.14) 

and 

0.6.15) 

The  characteristic  equation  along  the  particle-path  Mach  line  in 
the  solid  phene  la  either 
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Time  Period  1: 


K?„(dT»)n■,  Kt,  (dT)umKZ, (^a  (3-6-16> 

or 

Time  Perlc  I 2: 

‘KZ^Vn  (3-6-17) 

The  coefficients  before  each  tern  ere  denoted  by  K.*,  the  super- 
script x is  the  number  of  the  characteristic  equation  and  y the  subscript 
designates  the  coefficient.  They  will  be  shown  in  Appendix  3. 

Ths  characteristic  equation*  (3,6.11)  - (3.6.17)  are  the  additional 
information  which  is  needed  to  relate  the  flow  properties  on  tbe  boundary 
of  interest.  They  ere  used  together  with  the  governing  equations  to 
solve  the  left  boundary.  Kora  liar. use  ion  concerning  the  usage  of  these 
characteristic  equations  will  be  fiver*  in  tne  numerical  solution  sections 
later. 
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3.7  Initial  «nd  Boundary  Condition* 


The  ays ten  of  governing  aquations,  after  tha  above  aimplifications, 
becomes  * sat  of  six  first  order,  coupled,  non- linear,  inhomogeneous 
partial  differential  aquations.  It  is  found  in  tha  study  that  tha 
eigenvalues  of  the  system  are  six  distinct  real  numbers  under  noo- 
fluldized  conditions,  (Sea  Appendix  3).  The  characteristic  directions 
are  described  by: 

Us  + Gj  , (■^f‘ui~cs  > J (3.7.D 

(ofi TUp+c  > (-Tct Uf " C 1 Up  (3-7-2) 

The  subscripts  X,  II  and  III  represent  the  right-running,  left-running , 
and  gaseous-path  characteristic  lines  in  the  gas  phase.  The  subscripts 
IV,  ? end  TT  represent  the  right-naming , left-running  end  particle- 
path  characteristic  Haas  in  the  solid  phase. 

Wbm  the  granular  bad  becomes  fluidised,  the  speed  of  sound  is 
s»ec  trcx#3ltt*d  through  the  dispersed  particles;  therefore,  c is 
equal  to  aero.  However,  tha.  speed  of  sound  in  tha  gas  phase  is  never 
ser.>;  heacGc  tq  (3.7.1)  remains  the  same  but  Eq.  (3.7.2)  changes  to: 

r*  (3-7-3) 

whfea  i > d • l>ue  to  this  redaction  of  six  characteristic  linos  into 
c 

four  characteristic  lines  when  the  granular  bad  changes  froa  non- 
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fluidized  to  f Jl  uidizod  condition*  ? tha  governios  aquation*  cx&  uot 
totally  hyparbrlic.  Th*  solution  of  thsa*  equttio-i*  bateau  quite 
complicated.  Th*  ooundary  conditions  nsuyt  be  wall  specified. 

Tha  total  number  of  boundary  conditio*?*  required  depends  upon 
th*  flow  condition!  at  tha  ignitor  end  of  tha  grannie-  bad.  (A«  fsr 
aa  tha  currant  calcuiatieu  ia  concern*!,  ch«  ccaafoer  ia  sealid  by 
tha  shear  disc  at  tha  down*  ".ream  and;  th*  ph*roevux>a  after  tha  rapture 
of  the  shear  diac  aro  beyond  tha  acop*  of  this  study.)  Ta  specify 
cot.dit1.ona  at  th*  igniter  and,  * separate  control  voluma  (Fig.  3) 
is  considered,  which  ia  located  r mediately  at  tha  an trance  of  the 
granular  bei.  Whan  this  entrance  central  volume  contain*  both  gat 
end  particles,  six  ordinary  differential  equation*  ara  used  together 
with  the  necessary  co^atibiiity  relatione  to  deteoaina  tha  average 
flow  propartiea  in  tha  eu  trance  volume  and  also  the  properties  at 
th*  boundary  between  th*  entrance  volume  and  the  granular  bad.  Th* 
compatibility  relations  are  derived  directly  by  transforming  the 
governing  equations  into  their  characteristic  feme  (Appendix  3) . 

When  th*  entrance  eactlou  contain*  gas  only,  three  ordinary  differential 
equations  describing  the  time  rate  of  change  c t flow  properties  cf  the 
gas  phase  are  solved  together  with  the  appropriate  compatibility 
relations  for  the  boundary  paroaeters  reeded . The  total  number  of 
compatibility  relatione  used  in  th*  calculation  depends  upon  the 
flow  directions  cf  tha  gtues  tmd  particles  at  eha  borud-ry  and  also 
the  fluidia&tirra  conditions.  The  boundary  condition*  at  the  shear 
diac  end  are  relatively  simple;  before  the  shear  disc  la  raptured, 
tha  gas  and  pas ticle  velocities  are  all  aero. 
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The  initial  condition#  required  to  solve  the  eyetes  of  equations 

m%:*  tha  initial  distributions  of  gaa  temperature,  praaaure,  velocity, 

particle  velocity,  propellent  surface  temperature  and  fractional 

porosity.  The  initial  distributions  of  the  above  parsaetere  can  be 
non-unitors. 


5? 


3.8  Empirical  Correlations 


In  addition  to  th«  sbjvo  governing  aquations,  the  equations  of 
•tat*  for  gas-  and  partlcle-*phas«s  must  b*  specif lad.  The  co-;  ’luma 
effect  becomes  important  »t  high  pressures,  so  the  Nobie-Abel  dense 
gas  law  also  ccu.lid  the  Clausius  Equation  (75) , 

P(“J 5""bJ  —RT  (3.8.1) 


is  used  sc  the  equation  of  state  for  the  gas  phase  The  stateaitut  of 
a constant  density  for  th-*  solid-propellent  particles  serves  aj  the 
equation  of  state  for  tha  particles. 

To  complete  the  theoretical  model  it  is  osceaaar/  to  specify 
several  empirical  correlations:  the  intergranular  strase  transmitted 

through  tha  pecked  granular  particles,  the  flow  resistance  due  fu  tha 
drag  force  between  particle  and  gas  phesas.  the  convactiva  hast 
transfar  coaffidant  and  the  regression  rates  of  the  solid  propellent 
particles. 

A Constitutive  Lag  for  the  Total  itreee 

Tha  constitutive  lav  used  for  the  total  stress  calculations  is: 


whars  i*  the  critical  porosity  above  which  **h**to  la  no  direct 
contact  between  particle#,  and  c is  the  speed  of  sound  transmitted  in 


granular  material. 


j 

i 


T actually  consiata  of  two  parta  (1)  F,  tha  hydroatatlc  pressure 
(11)  tha  normal  strata  dua  to  tha  intergranular  normal  atraaa . x» 
tha  aolid  phase , P la  not  daflnad  lnaida  tha  particle* , l.a.  P - 0 and 
If  partldaa  vara  liquid  droplata,  P la  daflnad. 

Tha  ralationahlpa  of  tha  apaad  of  sound  in  tha  partida  phase,  c , 
and  tha  atraaa  tanaor  ara  given  In  Appendix  6. 

B Drat  Correlation 

Frcn  tha  fluidization  condition  of  the  gas-particle  system 
In  tha  granular  bad,  tha  bad  can  ba  divided  into  a non- fluidized 
(packed)  region  and  a fluidised  (dispersed)  region. 

Due  to  the  Importance  of  the  drag  force  term  in  the  momentum 
equations,  an  extensive  survey  vas  conducted  to  select  tha  most  suit- 
able drag  correlation  for  tha  modal.  Er gun's  equation  (45)  has  bean 
widely  adopted  by  some  pioneer  investigations  such  as  Kuo(l),  Krler 

(5),  Gough  and  Zwarts  (3)  in  their  studies.  Since  Srgun'a  equation  la 

*• 

valid  only  for  1 < < *»°00  and  0.4  < $ < 0.65,  it  has  vary 
limited  application  to  the  combustion  of  granular  propellents  under 
high  convective  burning  sltuatlonr.  A series  of  cold-flow  resists® ^ : 


measurements  were  made  bjr  using  a cylindrical  chamber  of  the  sa 


inner  diameter  cc  the  eombustloo  chamber  and  the 


type  of  pro- 


pellent, but  under  non-fluldlsed , non-combustion  conditions.  The 
correlation  (7$  obtained  la 


(3.8.3) 


which  is  valid  for  1 < *■  24,000. 


For  the  fluidized  region,  the  expression  for  is  deduced  from 

Andersson's  expression  ($6)  which  is  valid  for  porosities  renging  from 

0.45  to  1.0  end  Re  from  0.003  to  2,000.  We  have 
P 


where  the  tortuosity  factor,  tr,  is 


Q45&04.Q94’ 

Q9 


(3.8.5) 


The  cross-section  factor  Z and  inertial  drag  coefficient  ere 
given  by 


z 


I 

Ztr(l-0)0'59 


(3.8.6) 


<3-8'7) 

Although  the  Reynolds  nuaber  range  is  not  vide  enough  to  co*'sr 
the  variation  in  the  overall  transient  process,  it  is  the  best 
correlation  available  la  the  literature  . Since  the  values  of  D 

V 

in  the  fluidised  region  are  substantially  smeller  than  those  of  ths 
nonfluidlzed  conditions,  any  inaccuracy  Introduced  by  extrapolating  the 
drag  correlation  for  a fluidised  bed  to  a higher  Reynolds  number  will 
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cot  significantly  influence  the  calculated  results.  The  dreg 
eorreletlon  for  the  nonfluidited  region,  however,  is  very  lnportant 
in  the  aoaentua  equations. 


C Heat  Transfer  Correlations 


For  convective  heat  transfer  calculations,  Denton's  fornula  £7) 
is  used  for  the  non-fluidlzod  region,  hence 


IHizM 


fPr* 


It  is  valid  for  le^  fro*  500  to  50,000  end  porosity  around  0.37. 
Although  the  porosity  range  is  very  narrow  for  this  correlation,  it 
represents  the  porosity  of  aoet  of  the  packed,  unburned  portico  of 
the  granular  bed  In  ccnbuetloo  exper taints.  For  fluidised  regions, 
hc  is  obtained  free  JLove-Claacton's  correlation  £#), 


(3.8.9) 


where  the  exponent  of  the  leynolde 


is  expressed  as 


(3.3.10) 


The  Hewe-daacton  correlation  is  val4d  for  10  < Is.  < 10  . 

P 

las ides  the  oonv active  beat  transfer,  the  radiative  heat  con- 
duction is  rep rerouted  by 


h,«s  -€*  crCnt)  t r,* 


(3.1.11) 


« 


where  e ie  the  emits ini ty  factor  i»  taken  aa  unity  and  c la  the 
P 

Stefen-Bolr.xaan  constant.  Therefore,  the  total  heat  tranafer  is 
given  by 


h tft)-he(t)  +hr*lft) 


(3.8.11) 


and  is  given  previously  in  the  long  form  by  equation  (3.3.13) 

ht  (t)«  he  (thtp6  [T(tJ  +Tps(t)][T*r+)+Tp5*ft)]  (3.8.13) 


D Burning  Rata  Lae 

The  last  empirical  correlation  needed  to  close  this  problem 
is  the  burning  rate  calculation.  Since  the  relative  velocity  between 
gas  and  solid  phase  is  very  liigh  in  some  regions,  the  erosive  burning 
effect  oust  be  considered.  The  widely  used  Lenoir-RobiJ lard  (77) 
semi-empirical  burning  rate  lav  is  adopted  here 

r,  - o P%K*h,  «•«•«> 

where  fe,,  the  local  saro-b lowing  convective  obtained  from  either 
Kq.  (3.8.8)  or  Eq.  (3.8.9),  depending  upon  the  fluidised  situation. 
The  erosive-burning  const  ait,  K#,  and  the  erosive-burning  exponent, 

8,  used  in  this  investigation  are  listed  in  Table  1. 


TABLE  1 


DATA  PS ED  IN  THE  KOMEBICAL  COMPUTATIONS* 


°p 

1.6  g»/ca* 

a 

P 

0-945  x 10“ * ca2/aac 

k 

P 

5.30  x 10”"  cal/ar-see- 

k 

0.265  x 10  1 ctl/c»-sac 

V 

0.44615  x 10  * sa/ea-*ec 

b 

1.26  C3SS/gB 

T(500#X) 

1.4 

T<3000*K) 

1.26 

Pr 

0.7 

\ 

0.41275  x 10'1  an 

t, 

w 

0.3995 

A 

0.4745  cm* 

•a. 

15.24  cb 

T. 

293.91*1 

% 

1033.23  g/cm* 

♦. 

0.3995 ' . 

0.1  e»#-*X/cal 

i 

105. 

S 

1.0 

• 

0.65 

*X»t*  Car  ywptUt  properties  tn  obtained  fna  Ballistic  March 
Laboratories  of  B.i.  Army. 


Tabl  - S60*K 

lign  - 630*K 

T.  - 2500#K 

hg 

Tf  - 'S31*K 

C.(r  ) - 0.1969 

a p 
ro 

r^(undsc«rr«d)  - 0.032943  ca 

a - 0.362  x 10'* 

• o 

(cn/aac)/(g/ca*) 

a • 0.2814  x 10"“ 

so 

(ca/*«c) / (g/cn2) 


cd  - 0. 


00 

cd 

- 0.30 

10 

n 

00 

• 0.8867 

n 

so 

- 0.8020 

a and  n arw  linearly  lncar>olacad 

Cron  a » a and  n , n 
oo  at  i»  ii 

according  to  tha  xaluao  of  C,  at 

tba  laaenacanaoua  surface  of  eha 

grain. 


3.5  Ignition  Crifrion 


In  the  present  stud",  tbm  mathemati., cl  forBul*  cions  governing 
the  flcv  phanooana  in  c'  a gr  Hauler  bed  have  already  been  e complicated 
problem.  It  will  be  hi  1; i\y.  to  have  « simple  ignition  criterion  in 
completing  the  mathematical  modal. 

In  tbie  study,  two  crit  Leal  temperatures  vers  specified  at  the 
spherical  pellet  surf acme  to  achieve  the  ignition  criterion.  The 
lover  critical  value  corresponds  to  ths  ablation  temperature,  while 
the  higher  one  represent i the  ignition  temperature  of  the  propellent. 
The  burning  rate  starts  from  saro  as  soon  as  particle  temperature 
reaches  the  ablation  tempt. ature  end  increases  rapidly  to  the  burning 
zees  value  when  full  ignition  is  reached.  The  convective  heat 
transfer  lrops  to  sero  in  a linear  manner  as  ths  particle  temperature 
Inciuoau  from  the  ablation  to  the  Ignition  temperature;  Theoretically 
these  two  critical  temperatures  should  be  dependent  on  pressure,  but 
this  relationship  requires  further  -tudy.  Thus,  for  the  purpose  of 
this  stutfy  (since  pressure  dependence  at  high  pz  assure  la  smell) , 
they  could  be  taker  as  constants  In  this  study. 


CHAPTU  IV 


NCMX&XC&.  SOLUTION  TZCHNXQU3S 

The  previous  chapter  boa  provided  the  physical  background,  basic 
assumptions,  governing  exjunctions  and  sow*  associated  empirical 
corralationa  for  thia  two-phase , reacting  flow  nodal.  In  this  chapter 
wa  address  eha  aacbod  to  determine  tha  solution*  of  tha  governing 
aquations.  Tha  coup  laxity  of  tha  ays tan  of  governing  equation*  make 
it  impossible  to  obtain  analytical  aolutiona  for  tha  system.  As  a 
result,  va  turn  to  nuwrlcal  method*  fer  chair  solution*. 

Altar  a owe  careful  considerations , the  explicit  finite-difference 
achaw  is  aaployad  in  the  interior.  Method  of  characteristics  to- 
gether with  sow  compatibility  relations  are  us  ad  to  solve  tha  boundary 
points.  To  as sura  tha  convergence  of  the  numerical  solutions,  several 
convergence  testa  have  been  made  and  proved  to  be  very  satisfactory. 

Due  to  the  steep  gradients  in  sow  flow  properties  resulting  from  the 
rapid  transient  combustion  of  the  propellants,  numerical  instability 
is  encountered  when  the  straight,  forward  two-step  ticheeyer  cchaw  is 
adopted.  The  modified  version  of  the  two- step  Richsmyar  schwa  with 
a swell  amount  of  explicit  artificial  viscosity  are  introduced  into 
tha  computations.  Several  stability  criteria,  such  as  ths  Couraac- 
Trledrlcka-Lmvy  (C.F.L.)  conditions (69),  Energy  Method  (7Qf.  the  von 
New nn  stability  criteria (7j1  end  tha  Keuristlc  Stability  Theory  (6i)  are 
applied  on  the  finite-difference  equations.  In  dosing  this  chapter 
the  computation  procedures  and  the  flexibility  of  the  computer 
progr—  ere  die  cussed. 
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4.1  Selection  of  Nucarlcal  Kathode 


In  .-rd«f  co  selva  chw  sst  of  govaming  aquations  given  in 
Chapter  III,  several  different  approaches  asy  be  taken.  They  are 
the  eat  hod  of  characteristics  ( 74),  the  various  implicit  (70)  and  explicit 
f irate-difference  method  (71),  analog  sidle  cion (79),  and  sons  others. 

The  net hod  of  characteristics  is  supposedly  the  noet  suitable 
method  to  solve  s system  of  hyperbolic  equations.  Kitchen* (ll'1  has 
used  this  Betood  (with  *oas  success)  to  solve  the  net  of  governing 
equations  Kuo  (1) developed.  This  approach  has  been  known  for  its 
extreme  accuracy  but  Involves  complicating  and  lengthy  progressing 
logic.  An  undue  number  of  iterations  are  required  to  locate  the 
Intersections  of  the  characteristic  lines  and  on  enormous  mount  of 
interpolations  and  extrapolations  are  involved  in  obtaining  the 
flow  properties  at  desirable  locations.  All  tfci computation* 
are  rime  eonsvsslng  sad  leading  to  in*neut*cies.  Du*  to  the  above  reasons 
the  method  of  characteristics  was  adopted  At  the  boundaries  only. 

Implicit  methods  have  no  restriction  no  the  domain  of  dependence 
is  the  time  specs  field.  This  method  allows  large  time  stops  in 
computation,  out,  since  tise  flow  properties  change  too  drastically 
between  each  mode,  this  method  was  mot  applicable  to  our  studies. 

Stsmetor  (®0  used  Implicit  mthods  in  treating  setae  unsteady  flow 
calculations  hut  Abbott  (®D  found  that  these  methods  violated  the  "law 
of  forbidden  signals"  inberwat  to  system  of  hyperbolic  equations. 

Mam  wo  have  mr rowed  down  the  field  of  possible  techniques  to 
that  of  sa  explicit  ethane.  Its  slap  11  el ty  and  applicability  fur 


digital  computer*  niuera  Its  use  advantageous.  However , on* 
inherent  disadvantage  Is  its  instability.  One  may  overcame  this  by 
special  treatment  of  the  non-linear  and  inhomogeneous  terms. 

In  this  technique,  the  partial  derivatives  of  the  dupecdent 
vsrlables  are  approximated  by  difference  equation.  The  system  of 
partial  differential  equations  Ls  thus  transformed  into  a sal  of 
algebraic  equations . Different  investigators  use  different  methods 
to  tepresen’;  the  partial  derivatives,  thereby  arriving  at  different 
numerical  results.  Iu  order  to  choose  which  finite-differencing 
schema  tc  follow,  one  has  to  notice  two  impott&at  criteria,  namely, 
the  convergence  and  the  stability  of  the  method. 

Besides  hundreds  of  reports  and  papers,  detailed  discus a ion  cn 
finite  difference  methods  is  contained  in  textbooks  by  Forsythe  end 

'82).  Snfph(83l  Ragche  (71).  Vlchn*^*tskv  i,74).  and  Sichtmver  & Morton 
<70)L  Moretti (so,  84,  85,86)  has  also  done  extensive  work  on  this  topic. 

He  rakes  a clear  distinction  in  the  treatment  of  boundary  and  interior 
points  and  suggests  that  th*  authod  of  characteristics  be  used  at  the 
boundary  (exclusive  of  the  method  employed  for  the  interior)  . 

Instability  may  be  generated  if  the  boundary  has  not  keen  properly 
treated  and  will  cause  problems  in  the  interior  coo. 

Numerous  investigators  have  done  si&ilar  raaearch  to  this  present 
study,  they  are  East  and  McClur-&(4);  Fisher  end  Graves^);  Gough  and 
2warts(3)?  Gough  0-3?  Frier,  Van  Tasaell,  Pajan  and  VarShaw(5) ; 

Kitchen#  and  Kuo,  Vichnevetsfcy  and  Suaaserf  i#ld(8) . They  all  used 
a slightly  different  numerical  scheme  to  solve  their  governing  equations. 

Based  on  Che  valuable  axperlenco  of  reaeareners  in  the  field  of  granular 
prop?' lent  combination  we  have  mads  our.  choice  of  numerical  t shames. 
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The  method  of  characteristics  together  with  cor*  compatibility  relations 
have  been  adopted  in  the  bounded  calculation* . As  explicit  method  for 
the  interior  ha*  baan  u*ed,  namelv  the  modified  two-step  Mchrmyer  scheme 
incorporated  with  predictor-corrector  calculations.  In  using  the  pre- 
dictor-corrector technique,  the  non-linear  terms  are  treated  in  an  explicit 

t 

manner.  This  combination  of  the  method  of  characteristic*  and  finite 
difference  equations  has  bean  shown  to  ba  vary  efficient  in  the  present 
system  of  partial  differential  equations. 


4.2  Tba  Slaita-Eifferanca  dcuMilon*  for  tha  Interior  Points 

for  interior  point*  tns  modified  two-snsp  Richtoyer  scheme  wos 
used  with  predictor-corrector  calculation* » By  nh*  cwo-stup  modified 
Rich toy tr  schema  we  scan  a scheme  1a  which  tk*  diffusing  scheme  (8  7} 
and  the  modified  Leap-frog  scheme  6 2)  are  used  for  alternate  ulae 
cycles  In  the  following  manner. 

In  general,  the  system  of  differential  equations  are  represented 

by 


da 

ax 


(4.2.1) 


where  u Is  a vector  whose  components , m In  number , ere  'the  dependant 
variables  and  A In  this  case  Is  an  a by  m matrix.  To  simplify  the 
notation,  A is  consider sd  to  be  a constant  matrix.  The  system  is 
hyperbolic  if  k has  all  real  eigenvalues  and  m linearly  independent 
eigenvectors.  In  this  study  A was  found  to  consist  of  all  real 
eigenvalues,,  so  we  concluded  that  our  system  of  aquations  was  of 
hyperbolic  type  I represents  the  inhomogeneous  taros  appearing 
in  the  system  of  equations.  The  finite  difference  diffusing  scheme 
using  central  differences  in  specs  and  forward  differences  in  time 
is  employed  for  the  first-step  calculations.  The  difference  algorithm 
for  the  predictor  and  corrector  calculations  la  given  below. 

The  following  notation  was  used:  if  u (x,  t)  ia  any  function  of 

one  apace  variable  end  time,  denotes  the  finito-difference  approxi- 
mation to  u (i  A x,  j d t) , where  Ax  and  At  are  the  dimension*  of  the 
unit  EswPh  in  the  x,  t plane  and  i and  J have  Integer  values. 

The  first '«eep  diffusing  scheme  becomes . 
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Predictor: 


3U  _ 0.5(UJ+, 

at  “ 

.&L  <m  Ik^iL"  y^r.'  1. 

ax 


(4.2.2) 


Corrector : 


3u  tuf-nsbJi^uL)] 

at  m 


ax 


l (U.2.5) 


Th«  non- linux  coefficients  and  ichooogeneciis  seruat  v«re  tree tod 
explicitly  In  the  predic  tor-cox  cer.tor  technique. 

Using  Eqs.  (4.2.2)  end  (4.2.3),  Eq.  (4.2.1)  is  eKpr.osaed  In  another 
fora  as  follows, 

uf 1 - y (U* < * U w ) " ^ (***»  " (4.2.4) 


+ AtI 


i 

i 


•CfT  A[U-9)[iLUr  UU) 

+fle«J^uw-i>At[(i-e) if  j 


(4.2.5) 


Sqs.  (4.2.4)  and  (4.2;5)  srs  explicitly  Ae  fired , s*nnn  each 
aquation  determines  on*  unknown  in  tarns  of  old  quantities  uut. 
u£  represents  th«  predicted  quantity  iron  the  lstact  itaratius* 
of  tha  pr*di“.tor  calculations.  Whan  I • 0,  Sq.  (4.2.4)  reduces  to  th* 
original  diffusing  echos*  (which  has  second  or  dor  acen*ae>).  Tha 
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modified  Leap-frog  scheme  us lug  central  differences  in  both  space  and 
tlx*  is  employed  for  Che  feecoud-step  calculations.  The  Crsok-Nicolaon 
perimeter  (73)  8,  which  lies  between  0 end  1,  is  used  to  help  stabilize 
the  nvsmerical  solution.  The  difference  algorithm  for  Che  second-step 
in  the  predictor  and  corrector  calculations  is  given  belov. 

Predictor; 


3u  (at1  -ui) 

3t  2 At 

J+l  ]♦!  \ 

3U  -UU) 


ex 


z 


(4.2.6) 


Corrector: 


3u 

mrnmsmmm  Ml  ■ 

3-1  * .# 

L 


3 
o w 


au  r n-e)(uL-uL)+&(uJ**-u£*‘ 

ax  4ax 


//"'is 


(4.2.7) 


Using  EqSv  (4.2.6)  esc’  (4.2.7),  Eq.  (4.2.1)  can  be  expressed  in  the 
following  fore, 

-u  f **-  u!-^a  (of,  - u )+24t  if 


uf*  - ui a [a-eHul.-uL) 


*e(ut:?-utt‘)*(utt -«SS>J 
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Again,  Eq«.  (4.2.8)  and  (4.2.9)  ara  expressed  explicitly.  A 
•at  of  conaarvatlon  of  aaaa  aquation*  is  tha  gas  phasa  ara  iaplaaented 
is  tha  finite-different*  fora  using  tha  abova  aantlonad  nuaarical 
schema  and  ara  shown  is  Appendix  3. 

Tha  abova  predictor-corrector  tachsique  has  baas  us ad  to  overcome 
tha  nonlinearity  of  thaae  partial  derivatives  sad  tha  ishoaogaoaous 
ran*.  Is  tha  pradlctor  calculations,  tha  pravious  tins  quantitias 
ara  us  ad  to  avaluata  tha  coefficients  bf  tha  spatial  darlvatas  and 
inhomo par sous  tarns.  Tha  solutions  datarminad  from  thasa  calculations 
ara  avsragod,  according  to  8,  with  tha  pravious  tins  quantitias. 

Thasa  avaragad  valuas  ara  usad  in  calculating  the  non-linaar  co- 
affidants  and  inhonoganaous  tans  in  tha  corractor  calculations. 

This  convsaelonal  pradictor-corractor  tachniqua,  caused  nuaarical 
instabilities  in  both  tha  interior  nodes  and  tha  boundaries.  As  a 
result,  a snail  anount  of  artificial  viscosity  (8$  was  introduced  in 
the  calculation  to  filtar  out  tha  high  frequency  eoapononts.  This 
was  dona  only  in  tha  2 second-step  calculations,  sine*  ths  first- 
step  calculations  are  rather  stable. 

Following  tha  isp lamentation  of  this  two  atap  nuaarical  schaaa 
with  pradictor-corractor  techniques  and  an  artificial  viscosity  in 
the  second  stop,  a significant  iaprovaaant  of  atahlllty  has  bean 
achieved;  tbs  predictors  and  corractors  art  vary  dose  nuaaricslly. 
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4.3  Nu—rical  Treatment  of  ths  Boundary  Condition* 


The  boundary  points  art  calculatad  by  using  physical  and  corn- 
patabllity  ralations  to  achitvt  a high  dagrtt  of  numerical  stability. 

Tha  compatibility  ralations  us ad  can  also  ba  ragardad  as  ths  axtranaous 
boundary  conditions(l)  which  ara  naadad  whan  tha  cantral  dlffaranca 
tachnlqua  is  usad  for  spaca  darlvatlves  in  tha  hyparbollc  partial 
dlffarantlal  aquations.  In  this  study,  thasa  additional  boundary 
conditions  ara  a sat  of  charactaristlc  aquations  darivad  aarllar  in 
Chap tar  III. 

Tha  total  nuabar  of  boundary  conditions  on  the  left  and  of  tha 
chamber  depends  on  tha  nuabar  of  physical  boundary  conditions  plus 
tha  axtranaous  boundary  conditions  available.  Tha  nuabar  of  conpatlblllty 
characteristics  aquations  usad  depends  on  tha  flow  directions  of  tha  gas 
aaj  particles  — ths  fluidization  conditions  at  ths  boundary.  Eight 
sets  of  lsfi  boundary  conditions  havs  bean  considered  in  this  study. 

Tha  right  boundary  conditions  ars  simple,  tha  gas  and  particls  velocities 
ara  aero  before  tha  shear  disc  braska.  It  is  necessary  to  detetmlna 
how  tha  characteristic  aquations  worksd  on  a t-x  diagram  before  shoving 
tha  dlffaranca  approximations  that  apply  along  tha  charactaristlc  11ns. 

When  both  u^  (x^,  t)  and  u^  (x^»  t)  srs  subsonic,  and  ars  flowing 

In  tha  positive  direction  and  ths  bad  is  aon-fluldisad,  (LBCJ,  ) l.a.  Left 

Boundary  Condition  - Casa  1 , tha  Mach  Unas  at  ths  boundaries  ars 

shown  In  Fig.  4.  In  this  cass  there  era  two  left-running  Mach  lines, 

II  and  V,  seat  out  from  A^  sad  A#  passing  through  point  B at  ths  left 

boundary.  Tha  II  charactaristlc  aquation  with  a slops  • u^  - c 

amd  tha  ▼ characteristic  aquation  with  a slope  (4§)_  » u - c. 

at  v p 

Similarly  thare  ara  four  characteristic  Unas  passing  through  point  B ’ 
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(«_.«*>  (x 


,t-)  t*,.**)  (V' 

Figure  4 Characteristic  birectiane  at  the  Two  Boundaries  of 
the  Granular  Bed  when  the  Cae  at  x.  Flows  free  the 
Entrance  legion  into  the  Fropellent  Bed 


C^.t4*1)  (x^.t**1)  <V  •r3n> 


r)  <«#.r)  (x^*3) 


Figure  5 Characteristic  MreetloM  at  the  Two  lento  \ 
the  GtmmIit  Bed  when  the  Gae  at  x.  Flea*  free  the 

^ a _ _ * 


at  cha  right  boundary.  Thay  ara,  charactarlatica  linas  I,  III,  IV 
and  VI,  which  hava  alopaa  (^fjj.  - ug  + cg,  (^f)in  - ug,  (^)ly  - up  + c 
and  (j|)vl  “ up  raapactiaaiy. 

Tha  nunbar  of  unknowns  In  Cha  antranca  ragion  ara  six,  thay  ara 
u , u . 4 , T , P arid  T Thaaa  unknowna  ara  datarainad  aaally  by 

*c  Pc  c c c P*c 

tha  govaraing  aquation,  darlvad  in  saction  3.5.  With  tha  abova  qumtltlaa 
datarainad,  vs'coablntd  than  with  tha  two  compatibility  ralatlona,  givan 
by  Eq.  (A-3.46)  and  (A-3.58)  and  aotaa  othar  assumptions,  tba  alx  un- 
knowna at  point  B',  u , u , <fr,  T,  P cad  T ara  datarainad. 

3 P P* 

Sinca  tha  antranca  ragion  ia  «o  doaa  to  tha  laft  boundary  points, 
tha  following  assumptions  vara  aada  to  facilltata  tha  uaa  of  tha  two 
compatibility  ralatlona  Eq.  (A-3.46)  and  Eq.  (A-3.58),  thay  ara 


4 i /v 

Ai> 


a l 

W 


J I 


(4.3.1) 


(4.3.2) 


P(Xi.,t)  Pe("t)  (4.3.3) 

Fraai  eharactarlatlc  aquation  V,  Iq.  (A-3.58),  with  tha  assumptions  givan 
by  Eq.  (4.3.1)  wu  obr aland  tha  poroalty,  4,  at  paint  1 by  dir  net  aub- 
atitutiaa. 


76 


The  compatibility  relation  given  by  characteristic  II  i*  integrated 

along  its  charecteriatic  direction  to  yield  one  additional  condition,  u^. 

This  equation  ia  solved  simultaneously  with  the  other  two  physical 

boundary  conditions  Eqs.  (4.3.1)  and  (4.3.3)  and  the  calculated  4 ct  B. 

The  firs:  two  conditions  further  involve  two  tore  unknowns,  tne  ch&aber 

pressure  J (t)  and  the  chamber  particle  velocity,  u (t).  The  equations 
c pc 

involviac?  (t)  and  u (t)  are  given  by  Iqs.  (3.5.15)  and  (3.5.6), 

"•  Pc 

together  with  Eqs.  (4.3.1),  (4.3.2)  and  (4.3.3)  form  a closed  system. 

The  characteristic  equation  XI  can  be  solved  simultaneously  with 
the  other  equations  by  direct  substitutions.  Numerical  oscillation* 

/irs  observed  when  the  coefficients  of  the  differential  terms  ara 


evaluated  at  A (Fig.  4).  In  order  to  overcome  the  numerical  oscillations 
a 

in  the  boundary  calculations,  ths  sophisticated  fourth  order  lm#e- 


latte  integration  method  (1X4)  was  employed  (1$  to  carry  out  ths  in- 
tegration along  ths  laft-rumlag  charset  eristics.  In  ordar  to  psrform 
ths  1X4,  characteristic  equation  II  is  reconstructed  to  form  an 
ordinary  differential  equation  in  tb*j  following  msmnctr. 


J tdtk+KZmAt]/( KJj)  (4  3 4) 


Dm  preeeure  difference,  particle  velocity  diff 


, pereelty  diff 


ernes,  at  t Be  right  fcmd  side  of  Sq.  (e. 3. 4)  ere  the  properties  at  1 
ths  linearly  JmtarpoiLated  quantities  at  A(.  After  ths  foerth-o'der 


tge-Katta  iategratltti  is 


, ths  gas  velocity  at  the  life 
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(T,  + ?.Tt+ZT5  + T+) 


(4.3.5) 


Khar*  tha  coafflcianta  T , T , T «od  T ara  calculated  froa 

raj  i 

-iW-  KJkWPi  -Ki,  !„.(<*&  -KJ  lBl(«P)»-KfMUitKK5W 
Tj-(-  K?  l-u  LjdvJi-Kjl^oPk'K  Ll^t^KljU) 

(4.3.6) 

fov  Predictor  itaoa ; 

T4.-eK:i^(dP4-k^u(dufi-K;uw)»v*KLi,^t)tK2,i«) 

for  Corrector  a caps  > 

TW-KjUdPj.-K^I.W^pk-KjUCd^KtU^tyifKJjIa) 

Touparatura  of  cba  (M , T,  la  found  by  using  tha  gaa  static 
taaparatura  at 


Onca  tha  boundary  values  at  B hare  boon  determined,  tha  now 
location  of  A*  at  tha  next  tiaa  level  la  caleulatad  from  tha  a lope  of 
tha  left-running  characteristic  lino  at  tha  currant  tiaa,  l.a. 


7b 


(4.3.8) 


f 


j 


Similarly,  the  location  of  is  given  by, 


j*i  j 

**■  "^A, 


4 * 

T W 


(4.3.9) 


Tbs  positions  of  A and  A vary  slightly  irith  r*9p*ct  to  tins,  so 
2 f 

the  locations  of  A and  A are  not  involved  in  the  iteration  procedures. 

2 S 

The  above  calculations  are  shewn  in  a fir*  chart  in  Tig.  6* 

The  treatment  of  the  right  boundary  conditions  is  simpler,  since 

no  right  control  voIusm  is  considered,  is  general  we  only  nave  sis 

unknowns  to  solve,  saraly  «f,  u?,  T,  6.  t and  T^.  Bex  ore  the  stain* 

less-steel  dlephregn  breaks , the  right  and  of  the  ccwbuetim  chenber  is 

cloned,  the  gas  velocity  u and  particle  velocity,  u , ere  always 

B r 

egual  to  seru,  these  conditions  are  given  by 

U9(X«,t)*0  (*.3.10) 


UfiXm,  t)“0  ' (*.3.u> 

The  other  fetor  unhnaana  am  found  by  the  four  eLaractarlscics , Z, 
IUV  IV  cad  I.  fines  the  gas  *Ooeity  is  very  «cll  in  the  neighborhood 
of  the  shear  dine,  tfcr.  eolation  la  quite  stable.  To  reduce  sfcs  Mount 
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of  computation*  „ BK4  lntagradon  uatbod  la  not  snployad  hare.  To 
obtain  tba  porosity,  at  B'  w»  diractly  substitute  tho  sore  partldo 
valodty  Into  duiractarlatlc  IV,  |q.  (A-3. 46).  Using  tbs  s«ro  gas 
velocity , ssro  partldo  velodty  and  tbs  calculator  porosity  at  B ' , 

«s  substltuts  into  characteristic  aquation  X,  Eq.  (A-3. 40) , tba  gas 
praasurs,  F at  B'  Is  found  <■asi.11  scaly.  Slallarly  tbs  gas  tenpnrature, 
T,  Is  ob flood  by  substituting  tbs  calculated  prsssurs  at  B ‘ and  otbsr 
known  quantltlss  Into  cbaractarlstlc  aquation  III,  Eq.  (A-3. 32). 

Fro*  characteristic  aquation  VI,  Eq.  (A-3. 65)  tha  partlclt  tsnpsraturs 
Is  obtalnad.  Tbs  above  treatment  Involves  only  dir  act  restitutions, 
so  tba  right  boundary  values  azs  a sally  obtalnad.  Tba  cbaractarlstlc 
aquations  I,  III,  TV  sad  VI  ara  racastad  as  foil  tars , 

Cbaractarlstlc  aquation  I baconaa 


- Kj  (°  K?m  (dt)^ 


(4.3.12) 


Cbaractarlstlc  aquation  III  baceaas 


K5,sti,  (4,3.13) 


Cbaractarlatle  equation  IV  baconsa 


Cbaractarlstlc  aquation  VI  be const 


(d  Tgg  T)* + Kg,  Cdt  ^ 


(4.3.13) 


U 


Eqs . (A-3.12)  - (A-3.15)  are  coded  directly  into  the  computer 
program  in  the  above  manner.  (See  Fig . 6a  ) 

The  lications  of  A , A , A and  A for  the  next  time  level  are 

1 S 4 S 

determined  from 


(4.3.17) 


(4.3.13) 


(4.3.19) 


Since  the  above  locaCioca  change  vary  elightly  with  respect.  to 
time,  they  are  calculated  in  advance  from  the  elopee  of  the  character* 
iwtic  lines  at  the  previous  time. 

When  the  gas  velocity  changes  sign  at  the  left  end,  the  particle 
velocity  still  remains  positive  and  the  bed  is  non-fluidized  (LBC5)(9Q). 
As  shown  In  Fig.  5,  there  ar*  three  characteristic  lines,  emanating  from 
A*,  A^  and  Ag  passing  through  the  left  boundary  point  2.  The  character- 
istic llxias  at  the  right  are  same  as  those  in  Fig.  4,  calculations  at 
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B'  remain  Che  same. 


As  in  Che  previous  case,  all  Che  unknown  quantities  at  Che  chamber 
arc  firsc  calculaced  from  Che  governing  equaCions.  Some  assumptions 
are  made  Co  facilitate  the  use  of  Che  three  compatibility  relations, 

Eqs.  (A-3.46),  (A-3.52)  and  (A-3.62).  They  are  Eqs.  (4.3.1)  - (4.3.3)  , 
Porosity,  , at  boundary  point  B1  is  obtained  from  characteristic 
V,  The  fourth V'/rder  Runge-Kutta  integration  method  is  applied  to  both 

r 

compatibility  cond  >.ior/,  aJ;ng  characteristic  III  and  V.  Gas  tempera- 
ture, T,  is  deteralnaa  fro  a characteristic  III  and  gas  velocity,  u is 

o 

from  characteristic  V.  After  all  the  unknowns  at  boundary  point  B' 

have  been  determined,  the  above  calculations  are  shown  in  ?lg*  7 . 

> 

tha  locations  of  A and  A are  determined  using  Eqs.  (4.3.8)  and  (4.3.9) 

2 S 

and  new  location  of  A^  is  given  by 


(4.J.2G) 


Besides  the  above  mentioned  LBC  (Left  Boundary  Conditions),  six 
other  left  boundary  conditions  havs  bean  considered  depending  on  the  sign 

of  u , u and  tha  fluidization  conditions,  this  is  shown  in  Appendix  7. 

S P 

Tha  extraneous  boundary  values  calculated  in  the  above  manner 
converge  vary  quickly.  Since  tb<(  boundary  conditions  have  been  well 
treated,  the  amount  of  error  propagated  into  tha  interior  grids  am 
small. 

The  combination  of  finita-diffarance  and  the  method  of  character- 
latlcs  has  been  shown  to  ba  «n  appropriate  method  to  handle  this 
type  of  hyperbolic  pevtial  differential  aquations. 


4.4  Convergence  Teats  and  Error  Analysis 


It  Is  a formldabla  task  for  the  formal  error  estimation  and 
stability  analysis  on  the  numerical  algorithm  selected  for  our  intricate 
system  of  coupled,  non-linear,  Inhomogeneous  partial  differential 
equations.  This  task  is  compounded  by  the  complicated  boundary  con- 
ditions. In  order  to  bypass  these  nearby  Insuperable  difficulties, 
various  convergence  tests  are  performed  to  verify  the  numerical 
solutions  converge  to  the  genuine  solutions  of  the  original  differen- 
tial equations. 

The  convergence  tests  are  performed  by  varying  the  mesh  sizes  In 
time  and  space  increments.  The  solutions  change  slightly  when  the 
mesh  sizes  for  time  and  space  ere  reduced  by  e factor  c:  four.  Within 
the  present  scope  of  this  two-phase  combustion  problem,  our  numerical 
solutions  are  practically  unaffected  by  changing  the  mesh  sizes. 

Dus  to  the  feet  that  convergence  of  the  numerical  solution  depends 
on  how  well  the  problem  is  posed,  various  tests  are  performed  on  the 
sensitivity  of  tha  solution  by  changing  the  initial  conditions, 
boundary  conditions  and  some  other  important  parameters  such  as  tha 
initial  propellent  particle  size,  the  fractional  porosity,  the  burning 
race  lew,  the  heat  transfer  correlation,  the  drag  correlation,  the  con- 
stitutive lav  of  intergranular  stress,  the  ignition  criteria,  the 
ignition  mess  flow  rata  etc.  In  ell  these  tests,  the  numerical  solution 
are  bounded  and  change  only  slightly  for  small  perturbations  on  the 
parameters.  The  existence  of  the  neighboring  solutions  is  proved  end 
the  convergence  of  the  numerical  solutions  is  therefore  assured. 


4.5  Stability  Criteria 


When  the  partial  differential  aquation*  ara  represented  by  tha 
finiti-dif farenca  aquations,  high  frequancy  componants  will  not  ba  traatad 
accurataly.  This  inaecurata  traataant  of  tha  high  fraquancy  componants 
laads  to  tha  numarical  stability  conaidarationa . All  fiulta-di.ffaranca 
schaaas  poasas  high  fraquancy  componants . If  thasa  coarponants  ara 
dampad  or  don't  grow  rapidly  with  tins,  tha  schaaa  is  stab la.  On  tha 
othar  hand,  if  thasa  componants  grow  rapidly  and  finally  dominata  tha 
calculation,  tha  achana  is  said  to  ba  unatabla. 

In  gcaaral,  tha  finit*-diffaraaca  sy stans  cannot  ba  aaaaaaad  for 
instability  of  non-llnaar  sys tarns.  Sonatinas , tha  stability  crltaria 
for  tha  non~lisear  caaa  is  studied  as  if  it  is  a linear  system. 

Several  different  methods  of  stability  analysis  ara  known  in  this 
field.  Ve  selected  tha  following  methods  to  evaluate  tha  stability  of 
our  numarical  schema,  they  era: 

A.  tha  Courant , Friedrichs  and  Levy  Conditions  (C.F.L.)  (69), 

% , 

1.  tha  Energy  Sathsd  (7Q), 

C.  the  vetv  Neumans  stability  analysis  (73), 

j).  Heuristic  Stability  Theory  of  'Hirt(6$. 

A.  The  stability  condition  for  tha  two~step  modified  Klcbtnyar 
sir  hens  when  T.  • 0 in  Eq.  (a. 2.1)  and  the  system  of  partial 
differential  equations  Is  hyperbolic  is  that  the  value  of  jj~  should 
not  ba  greeter- than  tha  largest  eigenvalue  of  X.  Thi«  i*  shown  by 

^fmax|u,ic,|]sl 
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where  is  given  by  Eq.  (A-4.13). 

It  should  b«  noted  chat  the  selection  of  At  is  very  critical 
for  this  stability  condition.  This  criteria  is  used  as  a 
basic  requirement  when  selecting  the  size  of  At  and  Ax  in 
the  numerical  calculations. 

B . Energy  Method 

This  method  is  closely  related  to  an  energy  conservation 
principle  of  the  differential  equation.  In  many  cases  the 
conserved  quantity  Is  not  necessarily  related  to  physical 
•a*rgy,  but  generally  the  square  of  the  independent  variable 
as  shown  in  reference  (70), 

Two  subroutines  are  built  in  the  program,  mainly  to  check 
the  balancing  of  the  conservation  equations  in  both  the  gas 
and  particle  phase.  So  the  corrector  step  for  the  diffusing 
scheme  and  the  correct!  ? step  for  the  Leap  Frog  scheme 
solutions  are  checked.  The  results  ere  highly  setisfactory, 
aince  the  belance  of  ell  five  conservation  equations  are 
checked  closely. 

C.  the  von  Keumann  Stability  Analysis 

In  this  method,  e finite  Fourier  series  expansion  of  the 
solution  to  the  model  equations  is  made.  The  decay  or 
amplification  of  each  mode  is  considered  separately  to  deter- 
mine the  stability  condition.  Tile  amplification  factor, G, 
of  the  system  of  equation  lies  very  nearly  in  tha  unit  circl  e, 
‘.'he  mathematical  step  for  this  seal* ill ty  Analysis  is  shown 
in  .efermnce(7j). 
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1).  Hire's  Stability  Analysis 

A datallad  study  on  tha  numerical  aches*  usad  in  this  two- 
phase  combustion  problem  has  been  performed . Tha  finite- 
difference  aquations  using  our  prasant  numerical  schema  was 
applied  to  tha  conservation  aquations.  Each  term  of  tha 
dlffaranca  aquations  Is  expanded  in  a Tayor  aarlas.  The 
lowest-order  terms  in  tha  cxpaxuilon  must  represent  the 
original  conservation  equation.  All  tha  high-order  terns 
are  called  truncation  errors  which  account  for  the  compu- 
tational instabilities  of  the  finite  difference  method. 

All  the  diffusion  terms  era  collected  from  tha  truncation 
error  terms.  After  following  all  tha  numerical  steps  wa 
found  that  tha  affective  diffusion  coefficients  era  positive 
for  the  numerical  schema  wa  used,  thus  tha  computational 
stability  is  assured.  One  typical  way  to  obtain  the  effective 
diffusion  coefficient  for  the  numerical  solutions  of  a partial 
differential  equation  can  be  found  in  Reference  (66). 

Prom  the  above  four  stability  studies  we  have  concluded  our  numerical 
schema  it  very  stable. 


@f 


4.6  Artificial  Viscosity 


Thu  modified  two- step  Richtnyar  s'hema  la  believed  to  be  scab la 
when  char*  ie  no  inhomogeneous  and  tha  ays tan  of  partial  differential 
aquae  Iona  la  hyperbolic.  Whan  inhomogeneous  tana  ara  lndudad,  in- 
atability  ia  encountered. A aiailar  condition  vaa  obaarvad  by  Gough  and 
Zvarta  ( 3\  Krlar  at  al  ( 5X  Fortunataly,  nathoda  ara  available  to  over- 
come this  problem.  Gough  and  Zvarta  uaad  a damping  procadura  basad  on 
tha  smoothing  thaory  of  Shuman  (93)  to  circumvant  this  problem.  Krlar 
at  al  used  an  artificial  viscosity  r>  areas  tar  (92)  which  provided  damping 
effect  on  tha  high  frequency  components. 

To  overcame  our  problem  an  artificial  viscosity  e , which  provides 
damping,  has  bean  added  to  our  a<&ama.  This  parameter  is  applied  only 
to  tha  Predictor  and  Corrector  for  the  two-step  Leap  Frog  in  tha  following 
manna r. 


ur*-(  i -tour + (.  cr  ♦ 6 u;: 


i+i 


(4.6.1) 


After  this  artificial  viscosity  has  bean  added,  the  numerical  solutions 
have  no  mors  stability  problem  introduced  by  tha  inhomogeneous  terms. 
A3,jE  weighting  factor  (e  * 0.035)  applied  to  tha  neighboring  mesh  points 
would  eliminate  tha  high  frequency  components  effectively. 
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4.7  Calculation  Procedure  sod  Flexibility  of  the  Program 

The  theoretical  model,  described  by  a sat  oi  governing  aquations 
togathar  vitb  tba  additional  relationships  art  shown  in  Chap tar  III. 

Tba  simplified  governing  aquations  in  saetion  3.4srs  codad  in  a Fortran 
language  for  an  ISM  370-168  Digital  Coaputsr. 

A program  consisting  of  a rain  routlna  and  fifty-thraa  subroutinas, 
has  baan  developed  to  solan  this  two-phaee  eoabustion  problem.  Tba  overall 
calculation  procedure  is  shown  in  the  flow  chart  in  Fig.  8 . The  program 
is  basically  divided  into  four  sections : 

(a)  Predictor  for  the  first-step  diffusing  schema 

(b)  Corrector  for  the  first-step  diffusing  a chans 

(c)  Predictor  for  the  second-step  Leap  Prog  scheme 

(d)  Corrector  for  the  second-step  Leap  Frog  schema 

TVs  Qaaitaaae  4 f4»eS  Ss.  *alM«1fSa  bWa  lm#e  Va*m<l  awee  amm  J4  ama  meamnil 
me  *<»*•»  aw  mnamnnm  am  aesm  *naa  vwnsanna/  Wwaew* aawMO  | v9wen« 

to  calculate  the  right  boundary  conditions,  than  to  calculate  the  interior 
pcinea.  This  sequence  is  repeated  for  the  predictor  and  corrector  of  the 
two  steps. 

This  program  has  been  dasignau  to  be  as  gsnsral  as  possible  with 
the  following  ape  rial  fsaturss: 

(a)  h set  of  boundary  conditions  srs  coded  into  s separate  sub- 
routine to  facilitate  changes. 

(b)  With  Mall  codifications  , this  progrsn  could  handle  moving 
boondavy  conditions. 

(c)  All  empirical  correlations  are  coded  separately  and  could  bo 
changed  readily. 

(d)  Two  special  subroutines  are  developed  to  check  the  solutions 
fqt  the  corrected  lenrel  of  both  diffusing  and  lea?  Frog  schemes. 
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Fig uwi  S Goaoral  Flow  Chnx t tot  tho  Solution  of 
flowoniing  IquntioM 
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Thus  danonstrating  the  convergence  of  tha  solutions  by  the 
"Energy  uathod  " (70). 

(a)  To  inaura  accuracy  »all  raal  nuabera  ara  given  in  double 
precision. 

(f)  Tha  sash  sixaa  of  Ax  and  Ac  could  ba  aaaily  changed  for  obtaining 
initiation  at.daaired  axial  locaticna  and  apacific  tiaaa. 

(g)  Tha  progran  is  -written  ao  that  tha  solution  procaas  could  ba 
a topped  and  restarted  at  ny  nine  in  tha  calculation. 

(h)  Coonon  blocks  ara  coapoaad  of  functional  groups  of  progran 
par  aan tars. 

(i)  Meaningful  nonane latur a haa  bean  aaaignad  to  the  variables  for 
a any  identification. 

Cj)  The  main  progran  contains  only  call  atatanants . 
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CHAPTER  V 


COMPARISON  OF  THEORETICAL  PREDICTIONS  WITH  EXPERIMENTAL  DATA 
5.1  Experimental  Setup 

Experiments  to  verify  the  theoretical  model  ere  carried  out  using 
the  cylindrical  chamber  shown  in  Fig.  1.  The  geometric  dimension*  of 
the  chamber  ere  given  in  the  usi  figure.  Tbs  spherical  propellent 
used  for  the  tests  is  type  WC-870,  with  sn  svmrsge  particle  disaster  of 
0.8255  oa.  Some  physical  parameters  and  chemical  data  of  the  spherical 
propellent  is  given  In  Table  2.  The  particles  are  packed  in  a chamber 
15.24  cm  long  with  0.777  cm  I.D.  to  form  a propellent  bed  with  an  average 
weight  of  6.95  grots.  The  cylindrical  chamber  has  a wall  thickness  of 
3.937  cm  which  is  much  thicker  than  the  downstream  stainless-steel  shear 
disc  (0.81  oa.).  This  0.81  mm.  burst  diaphragm  provides  a fixed  boundary  and 
also  serves  es  e safety  valve  for  the  chamber.  Consequently  , there  is 
no  permanent  deformation  to  the  combustion  chamber  after  each  teat.  Due  to 
this  fact,  the  sssm  apparatus  is  used  repeatedly  ft-:  s series  of  identical 
tests  when  the  loading  density  is  kept  constant. 

To  measure  pressure  and  flama  front  speed,  high  frequency  Minihat 
prtaaure  transducers  (93)  and  ionization  probes  (Dynasen  CA-1040)  are 
equally  spaced  along  the  chamber.  (The  pressure  transducers  are  statically 
calibrated  by  a dead  weight  tester  at  Ballistics  Research  Laboratories  of 
Aberdeen  Proving  Ground. ) They  are  placed  at  the  prescribed  locations  marked 
by  Gl,  G2,  G3,  G4  and  G5  in  Fig.  1.  In  ordur  to  maximize  chamber  strength, 
the  transducers  and  ionization  probes  are  placed  spirally  along  the  chamber. 

A gaseous  pyrogen  Ignition  system  is  used  to  ignite  the  propellant 
bad,  and  is  shown  schematically  in  Fig.  9 . This  type  of  Ignition  system 
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TABLE  2 


PHYSICAL.  COMPOSITIONAL  AND  THERMO-CHEMICAL  DATA 

— 6f  MTOTflgmmn* 


Granulation 

Max  Parti cla  Diameter  • 0.0965  cm 
Mja  Partial*  PI  a— tar  ■ 0.0666  cm 

Partial*  Shape  ■ Spherical 

Gravimetric  Density  ■ 0.960  gm/ce 

Z Nitroglycerin  • 10.0 

Z Nominal  Nitrogen 

Content  of  nitrocellulose  » 13.15 

X Da tar rant  Coating  » 5.20 

Heat  of  Explosion  - 870  cal/gm 

Flame  Temperature  “ 2831*K 


♦This  WC870  PROPELLENT  NAS  MANUFACTURED  BY  OLIN 
CORPORATION,  WINCHESTER  GROUP 


15 


Spark  Plugs 


Figure  9 Scheaaric  Drawing  of  Igniter  Fjed  Line  and  Ignition  Chaaber 


vas  ehoaan  over  tin*  conventional  impact  Ignition  primer  because  it  is 
capable  of  achieving  both  high  reproducibility  ea  veil  as  varying  Igniter 
strength  and  duration,  these  features  are  not  characteristic  of  pallet 
igniters  (9 The  length  of  the  igniter  chamber  4.763  a art  the  inner 
diameter  of  th*  chamber  is  0.777  ea. 

Igniter  seas  flow  rate  is  determined  by  the  choked  flow  equation 
using  Ignition  chamber  praaaur*  and  n .taperaturs , and  in  in  Fig.  10. 

The  Igniter  strength  and  reproducibility  are  controlled  by  regulating  the 
gaseous  mixture  pressure  and  the  fuel-oxidant  ratio.  The  renctont  com- 
bination of  gaseous  hydrogen  and  oxygen  wee  chosen  for  the  igniter  because 
of  its  vide  flamability  limits  and  because  it  is  easily  ignitad  by  sparks. 

The  gaseous  reactant*  art  vell-airad  in  the  ignition  chamber  and  the 
gsi  fmmd  lines  are  dosed  before  ignition  by  spark  plugs.  iVo  spark  plugs 
axe  used  to  insure  a successful  ignition  in  every  firing. 

The  ignition  and  combustion  chsmbart  ara  separated  by  a nuiri- 
perf orated  convergone  nozzle,  vhi -h  consists  of  seven  holor  of  1/16  inch 
IB.  For  graeter  flow  efficiency,  the  dive  tor  of  the  holes  it  expanded 
to  3/32  inch  at  the  interface  between  the  nozzle  and  the  igniter  chamber, 
forming  e parabolic  entrance.  A detailed  description  of  the  noxxls  is 
shewn  in  ?J.§ . 11.  Tbt  ec.ec ream  tide  af  thf*  nestis  is  covered  and  tightly 
sealed  by  a tape  to  prevent  unburnod  react  sot  gas  :*  froa  entering  the  pro- 
pell  vt  bad  before  ignition . The  design  eS  ito  mulci-perfurar.id  nozzle 
la  intend ad  specifically  to  provide  s radially  uniform  heating  of  the  granule: 
propellant  bed,  to  e*^iev»  a one^viMansionAl  flow  ennfition  ia  the  axial 
direction,  end  also  to  provide  c o*»-dirwnai«nr>  crsspc.ct_.cu  nf  bad 

The  transient  data  obtain  ad  from  prrsaur*  trat/sducer*  and  ionisation 
plae  era  resold  ad  in  the  I.M.  ttafie  of  a multi -eturoal  tape  recorder  at  a 
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10  Ignitor  Mm*  Flew  Bat*  Flowing  into  eh* 

ftcsa&Ur  i«d  «i  a Function  of  tin* 


Figure  11  Schematic  Drawing  of  the  Mult i-Porf orated  Convergent  Woizle 


speed  of  60  ips  In  order  to  provide  totally  time-correlated  results. 

The  analog  data  are  then  recorded  on  a transient  waveform  digitizer 
(Bionation  Waveform  Recorder  Model  1015)  which  obtains  100  daca  points 
per  millisecond  of  transient  data.  The  data  is  then  expanded  by  three  orders 
of  magnitude  in  time  and  plotted  on  a conventional  x-y  plotter  at  a rate 
of  oue  data  point  per  second. 
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5.2  Discussion  of  Results 


Ths  theoretical  results  ar«  obtained  by  using  the  mathsnatical 
model  prsssncsd  in  Cbapcsr  IH.  In  actual  computer  coding  tha  simplified 
aquations  (saa  saction  3.4)  arc  usad.  Tha  solution  method  is  shown  in 
Chaptar  IV. 

Tha  pradictad  valuas  of  pressure-tins  tracas  at  various  locations, 
such  as  G2,G3,  G4  and  G5  art  shown  in  Tig.  12.  As  in  a typical  sat  of 
experimental  data  given  in  Fig.  13,  tha  prassuriaation  process  at  tha  down- 
stream portion  increases  faster  and  eventually  overtakes  tha  upstream 
pressure  trace.  Further,  tha  data  obtained  from  pressure  transducers  show 
that  tha  rata  of  prassuriaation  invariably  increases  in  tha  downstream 
direction.  In  Fig.  13,  tha  pressure  transducer  Gauge  1 (Cl)  was  positioned 
in  tbs  igniter  chssber,  and  bacsise  of  a choked  heckflov  condition  at  the 
iwasla.  It  was  found  to  lag  behind  tha  rise  In  pressure  at  Gauge  2 location. 
Tha  four  other  gauges,  62  to  G5,  were  positioned  along  tha  propellant  bed 
with  tm  equal  spacing  of  3.T?3  cm  (Sea  Fig.  1). 

Tha  pressure  traces  measured  in  tha  granular  bad  indicats  that  an 
upstream  pauga  (such  as  G2)  senses  the  pressurization  much  sooner  than  the 
gauges  at  downstream  positions . however,  the  rate  of  pressurization  is 
significantly  higher  for  the  downs tremr.  gauges.  It  is  also  interesting 
to  note  that  the  pressure  peak  for  65  io  slightly  ohoad  of  G2.  Thlo  Is 
frscmaaa  tha  depressurization  affect  due  to  shear  disc  rupture,  o vs r cones  ths 
pressurization  due  to  propellent  gasification  at  a gives  station,  the 
pressure- tins  traces  start  to  decline.  The  pressure  at  ell  stations  drops 
heck  to  one  atmosphere,  ending  e transient  interval  of  about  5 nsec. 

The  predicted  preasvre-tlne  traces  correspond  to  locate  on*  62,  G3, 

64  and  65  are  individually  super inpoaoe  eh  sane  experimental  datii  for  six 
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separata  firing*  (tested  under  tha  exact  condition)  ara  shewn  on  Figs. 

14,  15,  16  and  1/  respectively . tha  shadad  raglon  rapraaa-ta  tha  dagraa 
oi!  swatr.ar  of  tha  axparlau.nr.al  data.  Tha  solid  llnaa  raprasant  ths  theoreti- 
cal pradictlons  for  tha  r ■ : in  ; t a# aura  at  tha  ccrvaspondlng  stations. 

Sine*  tha  theoretical  calculation  Is  limited  t tha  tits*  interval  prior  to 
tha  rupture  of  the  shear  disc,  no  comparison  is  made  for  tha  depressurization 
processes. 

Fig.  14  shows  a major  portion  of  the  pradletad  pressure  tree*  at  G2 
lias  in  between  these  fixing  data,  for  tha  early  part  of  the  transient 
interval,  tha  calculated  value  of  pressure  is  locar  then  the  experimental 
data.  This  Is  probably  due  to  tha  satirical  correlation  used  for  the  drag 
tan  in  tha  granular  bad  of  inert  apharss  (Eqa.  3.8.3  or  3.8.4)  and  may 
result  in  higher  resistance  than  the  real  test  condition.  In  Fig.  15  and 
16,  tha  comparison  between  experimental  and  theoretical  prassnre-tim*  truce* 
reveals  tha  same  trend.  That  is,  tha  caitclacad  pressure  at  G3  and  G4 
begin  to  rise  at  a lator  tins  than  tn*  experimental  data.  In  Fig.  17, 
tha  predicted  pressure  rises  almost  at  ths  same  clmtv.  of  tha  experimental 
results.  This  gauge,  G5  is  particularly  important;  for  its  location  near 
tha  shear  disc  allows  us  to  examine  ths  strength  of  tha  pressura  spike  clore 
to  the  shear  disc. 

It  is  shown  dsarly  from  the  above  comparison  th.-uv  the  predicted 
pressure- time  variations  at  G2,  G4 , G5  gauga  locations  not  only  have 
the  right  magnitude  but  also  have  rees<9nabj.e  slopes.  The  pressure  rise 
times  are  also  oot  for  frra  the  experimental  measurements . Although  the 
comparison  at  G3  la  lees  satisfactory  than  those  at  G2,  G4  and  G5,  the 
•lope  of  the  pressure-time  traces  at  G3  la  not  far  froc  chat  which  wee 
measured. 
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13  Cospsrlson  of  Theoretically  Frs41cte4  53  Pruiun  Trees 
with  the  Coopoelte  fracture  Trace*  of  Six  Ixperioentel  Firings 
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Figure  16  Costparlsee  of  Thaorotieally  Prodiet«4  G4‘ Pressure  * * »ca 
with  tho  CwpMlti  Pressure  Treeee  of  Six  Experimental  Firings 
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Figure  17  Comparison  of  The ore tic ally  Predicted  G5  Pressure  Trace 
With  the  Composite  Pressure  Traces  of  Six  Experimental  Firings- 


In  Fig.  18,  the  calculated  flam*  spreading  rat*  is  compared  with 
tut  data  from  thra*  separate  firing*  with  ionization  pin  meeaurements . 

(Three  other  firing*  of  th«  maa  series  consisting  of  tn*  six  run*  men- 
tioned above,  had  no  flam*  front  measurement*).  The  calculated  results 
again  check  quit*  cloaaly  with  the  experimental  data.  Both  the  theoretical 
result*  and  experimental  data  indicat*  a significant  increase  in  flare 
velocity  as  the  combustion  wave  propagates  In  the  shear  disc  direction. 

Fig.  19  show*  pressure  distributions  calculated  for  various  tines, 
with  the  locus  of  the  ignition  front  (dashed  curve)  superimposed. 

Beginning  with  the  trace  (t  • 1.26  msec),  the  pressure  distribution 
quickly  dsvelops  s continental  divide  and  the  gradient  on  its  right  side 
becoaes  steeper  with  increasing  tine  from  the  increased  gasification  rate 
in  the  combustion  sons.  The  pressure  gradient  grows  to  the  left  of  the 
peak  due  to  the  reverse  flow  of  gases  from  the  granular  bed  to  the  ignite? 
chamber.  Thus,  vten  the  pressure  at  the  head  end  of  the  granular  bed 
be cones  higher  than  that  in  the  igniter  chamber , the  igniter  chamber 
serves  as  a ness  and  energy  sink  to  the  combustion  chamber.  The  resulting 
re-verse  flow  of  gases  also  makes  the  peak  more  pronounced.  Further,  the 
reverse  gas  flow  eventually  causes  the  pressure  trues  of  t • 1.54  to  cross 
that  portion  of  the  prassnra  distribution  at  t » 1.42  msec  (this  can  also 
be  seen  from  Fig.  14).  Max  the  end  of  the  chamber  the  locus  of  the  Ignition 
front  curves  sharply  upward  (Fig.  19),  dua  to  ths  extremely  rapid  pressur- 
ization near  tha  sheer  disc. 

Fig.  20  shows  ths  temperature  distribution*  (et  times  corresponding 
i Fig.  19),  with  the  locos  of  tbs  ignition  front  shown  by  n deshed  curve. 

The  advancement  of  tbo  combustion  wave  In  the  granular  bed  is  readily 
seen,  larly  in  the  transient,  the  praetors  In  the  system  is  lew,  so  in 
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order  to  heat  the  granular  propellent*  to  ignition  condition,  the  gee 
temperature  at  the  ignition  front  mutt  be  high.  As  the  pressure  riaes, 
the  heat  transfer  from  tne  gas  to  the  unbumed  propellent  becomes  more 
efficient  and  the  ignition  front  locua  descends.  Finally,  near  the 
end  of  the  granular  bed,  where  tb«  gaa  velocity  is  small,  the  convective 
heat  transfer  effect  decreases  considerably  and,  theiafore,  the  gas 
temperature  at  the  ignition  front  again  riaes. 

Fig.  21  shown  the  particle  tes^  rature  distributions.  In  this 
figure,  the  locus  of  the  ignltlra  front  is  not  shown,  because  it  is  only 
a straight  line  across  the  temperature  rcale  with  T equal  to  T ... 

pi  iDl 

The  rate  of  flame  front  spreading  is  observed  to  be  more  pronounced  at 
the  later  stage  of  the  transient  event. 

Fig.  22  shows  the  corresponding  gas  velocity  distributions . Early 
in  the  transient , Igniter  gas  flows  into  the  granular  bed  and  the  gas 
velocity  la  greatest  near  the  noxzle  end.  Later,  with  combustion  occurring 
in  the  granular  bed,  a peak  develops  in  the  velocity  profile  driven  by  the 
pressure  gradient  near  the  ignition  front.  Still  later,  the  pressure 
gradient  in  the  granular  bed  near  the  noxzle  causes  reverse  flow  of  hot 
gases  into  the  igniter  chamber;  therefore,  the  gas  velocity  is  negative  in 
the  left  portion  of  the  combustion  chamber. 

In  Tig.  29  the  particle  velocity  distributions  are  shown  on  the  sane 
seels  to  emphasize  the  fact  that  particle  velocities  are  in  general  much 
lower  chan  the  gas  velocity.  The  motion  of  the  particles  are  generated  by 
the  pressure  gradients  developed  In  the  granular  bad.  Later  la  the  transient, 
a fraction  of  the  particles  move  in  the  reverse  direction  near  the  noszle 
end. 

The  poroalty  distributions  in  Fig.  24  show  sons  interesting  phenomena. 
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Initially  they  ara  quit*  uniform,  but  vh«:t  combustion  rtarco  at  the 
noicla  and,  grain  motion  generated  by  tha  prassur*  gradient  ra-pacl:  the 

granular  bad  to  reduce  the  poroaity  ahead  of  the  ignition  front.  Behind 
tha  front  it  increases  as  a result  of  burning  and  grain  motion  in  the 
forward  direction. 

Quantitatively , tha  calculated  and  experimental  results  era  in  good 
egraeaent.  This  verifiaa  that  tha  simplified  aquations  in  taction  3.4 
ara  truely  valid. 

The  shear  diac  that  seels  the  right  end  of  the  chamber  is  calibrated 
to  withstand  a maximum  pressure  of  773.4  Kg/ at2.  The  numerical  calculation 
stops  whan  tha  pressure  near  che  shear  disc  reaches  this  critical  value. 

Tha  calculated  results-,  show  that  the  diaphragm  breaks  at  1.638  msec. 
Similarly,  in  tha  experimental  process,  combustion  was  terminated  imsadiataly 
afier  tha  shear  disc  raptured,  due  to  depressurization  o * tha  chamber. 
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CHAP  TEX  VI 


CONCLUSIONS 

An  investigation  has  been  conducted  for  a tvo-phaas.  granular  pro- 
pellent coabo*tlon  problem.  Several  important  conclusion*  arc  suamarlssd 
in  tbs  felloe  lag. 

1.  A successful  theoretical  modal  has  been  developed  to  describe 
the  combustion  of  mobile  granular  propellant*.  The  calculated 
pres cure- tin*  traces  end  flame  propagation  rates  agree  closely 
with  the  experimental  results  aessured  by  high  frequency  pressure 
transducers  and  ionisation  probee. 

2.  A modified  two-step  XI  ch  coyer  explicit  scheme  is  rap  loved  is 
solving  the  system  of  hyperbolic  partial  differential  aquations. 
The  solution  shows  that  this  method  is  stable  and  fast  con- 
vergent whan  the  boundary  conditions  ara  solved  with  the  method 
of  characteristic* . 

3.  According  to  the  theoretical  calculations  and  experimental 
observations,  soma  special  features  of  the  combustion  procass 
la  granular  propellents  ara  notad. 

(a)  The  rate  of  preeeurisetion  Increases  in  tbs  downstream 
direction  and  the  pressure  peak  developed  in  the  granular 

bad  travels  downstream. 

. (b)  The  flame  spreading  rate  increaoaa  significantly  in  the 
downstream  direction. 

(c)  A stronger  igniter  ceuaes  faster  flame  spreading  and  also 
e higher  rate  of  p reason satlen» 

(d)  Ignitor  gas  flow  compacts  the  bod  and  causa*  the  propellent 
grains  to  msva  in  the  shear  disc  direction. 


U» 


(«)  A » tho  pressure  p««k  develop#  In  the  chamber,  the  g ca  and 
particle  respond  to  Chs  pressure  gradient  by  moving  from 
center  toward  both  ends  of  the  chamber  simultaneously. 

(f)  lor  tightly  packed  beds,  the  particle  velocity  is  sig- 
nificantly lower  than  the  gas  velocity  before  the  rupture 
of  shear  disc. 
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APPENDIX  1 


Conservation  Equation!  In  the  Eulerian  Porn 
The  governing  equation*  for  both  the  gas  and  particle  phases  in 
the  Eul arias  fora  can  ba  deduced  directly  froa  the  governing  equations 
In  their  divergence  fora  given  by  Eq.  (3.3.3)  - (3.3.8),  assuming  A is 
constant. 

The  mess  equation  of  the  gas  phase  in  the  Euler las  fora  is 


The  mess  equation  of  the  particle  phase  is 


^ (A-1.2) 

The  aomantua  equation  of  the  gas  phase  is 

^A»ript(^at)-s(AsPv^^li^l!a-)  (A-i.3) 

The  momentum  equation  of  the  particle  phase  is 


The  energy  equation  of  the  particle  phase  is 

s & 

-4.+a»{<w9  *$$-& -htrr-T»>+/Hi  Dw. 

-c.T.agrt  u-i.3, 


APP1HDU  2 


^najjtical  Solution  of  the  Solid-Fhaa a Beat  Equation 
Tha  heat  equation  for  imburnad  spherical  pellets  at  a fixed 
location  is  givon  by  Eq.  (3.3.9X  which  can  ha  ’■'■tnaforaed  into  the 
Fourier  haat  aquation  for  tha  slab  caaa  (9$  by  defining 


Tr-Tr-T6 


(A-2.1) 


Tha  heat  aquation  becoaes 


(A-2.2) 


Tha  initial  condition  baconaa 

t-o:  Tjto.n-O 

The. boundary  conditions  becoaa 

r-o  : $*£(t.0)-0 
r-r*: 

mUt) 

Wo  now  propose  a function  for  T^*  (t,r)  so 


(A-2.3) 


(A-2.4) 


(A- 2. 5) 


Tima  Period  1 : 


Before  tha  thermal  wava  has  penetrated  to  the 


center  of  tha  sphere,  l.a.  5 < r (Saa  Fig.  25). 

p. 

Thasa  four  constants  C , C , C and  C ara  to  ba  daterminad  by 

0 12  i 

four  conditions.  Threa  of  these  conditions  are. 


( A-2.7) 


t;  rt,  r,.  - J J-o 


(A-2.8) 


and 


<A-2*9) 


plus  tha  aaoo thing  condition: 

I$rt,  rw).o 


(A-2.10) 


which  tends  to  aaka  tha  teapavature  profile  go  smoothly  into  tha  initial 

temperature.  Wa  define  r . 5 r >5,  and  after  soaa  algebraic  aanipu- 

6 P. 

lations  and  substitution  C , C , C , C ara  as  follows: 

« i » i 
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(A-2,11) 


t 


1 M t. 


r _ 3 nr  rr~Z 


2r>*+r|  -3rfnr 


(A-2.12) 


(A-2.I3) 


<i <jU __ 

2-rV+r<r--3r>„*r> 


(A-2.14) 


Substituting  (A-2.11)  to  (A-2.14)  In  (A“2„6)  va  get 

T!ft  r) .rvA . fr-rf jJ  „ 

' P'>uj«  ~ ( ? r-_  3 * rv  r* 

v ■ rg  - 5 • r, 


(A-2.15) 


for  r > r. 


differentiating  Bq.  (A-2.15)  once  we  have. 


cr+rMc-jaiL 

ar  * “ (zr*+rtxn.-rff 


(A-2.16) 


differentiating  Eq.  (A-2.15)  twice  we  have. 


*/  _ ZnSi  r r'-rt  i 


(A-2.17) 


When  r ■ r T * (t,r  ) - T * (t)  eubctltutlng  this  into  (A-2.15) 
Pf»  P P P* 


Tn'6^~W^Ll 


(A-2.16) 
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Eq.  (A-2.18)  axpraaaaa  a relationship  between  the  pellet  surface 

taaparature  and  tha  penetration  diatanca.  As  a raault,  thara  ia  raally 

only  ona  unknown  f(t)  [either  T * Ct)  (O ) In  Eq.  (A-2.15)  and 

thla  will  ba  obtalnad  frou  tha  heat -balance  integral. 

Xnta grata  Eq.  (A-2.2)  with  raapoct  to  r fron  rfi  to  r , wa 

p o 

hava 


Tha  BBS  of  (A-2.19)  can  ba  expressed  aa 

(A-2.20) 

The  LBS  of  (A-2.19)  can  bo  expressed  aa 


d r(i>.-rrfath 
* *ti  rtrP,+  r»  J 


(A-2.21) 


Aftar  aubstltutlng  Bq.  (A-2.20)  and  Eq.  (A-2.21)  into  Eq.  (A-2.19) 
and  boob  raarranglng  of  tana  wa  eom  up  with  an  ordinary  differential 
aquation  by  using  this  integral  nathod, 

Q BkZil  iZjjlkzSS.  (A-2.22) 

d T 

la-MpMattg  thla  aquation  In  tana  of  "y*  » r,t  la  raplaead  by  6 
In  Bq.  (Jr2.ll)  aa  hava 
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Differentiating  Eq.  (A-2.23)  with  reject  to  t and  rearranging 


we  obtain 


.apkztfr  JaX-  4L.7 

3T  Tib/-  L dt  (3 rfo-s)  at  J 


(A-2.24) 


Substituting  Eq.  (A-2.24)  into  Eq.  (A-2.22)  and  using  the  definition 

of  2 (t)  and  the  definition  of  T * we  finally  have 

P 


SR- 


&£. 

5Tr< 


iS&IbcS!  i he  ? j 


where  6 can  be  obtained  fron  Eq.  (A-2.23)  with  5 expressed  explicitly. 


rw- 


-5».CTwft>T.7 


[Tpa(t)-T.]+  J2thdL>[T(t)-TM(t)] 


(A-2.26) 


for  0 < 6 < r 


Eq.  (A-2.23)  is  not  veil  posed,  because  when  S • 0,  it  will  blow 

up!  To  bypass  this  initial  singularity,  the  initial  value  of  6 can  be 

baaed  ee  the  initial  values  of  T..  T and  h„  which  are  slightly 

pe  t 

different  fron  undisturbed  values. 

When  6 has  reached  the  center  of  the  sphere,  the  tent  Mature 
distribution  can  be  obtained  fro*  Eq.  (A-2.15)  by  setting  r^  • 0 


(A-2.27) 


This  parabolic  profile  serves  as  Che  initial  profile  for  the  tenperature 

distribution  at  Tine  Period  2,  t represents  the  time  required  for  the 

P 

thermal  wave  to  penetrate  to  the  center  of  the  pellet  in  ignition. 

In  the  similar  way  the  heat  aquation  in  the  spherical  particles  for 
Time  Period  2 is  derived  (after  the  thermal  wave  has  penatratad  to  the 
center  of  the  pallet). 

Sam*  aquation  as  Tima  Period  1 1»  used, 


Initial  condition: 

Tp(t»,t)«Te+  j^[T  ftp)-T»rtP)J 


Boundary  condition  1: 


(A-2.30) 


Boundary  condition  2: 


djt 

dr 


•ft.  r>.)  t)-Trt,rP,)j 


(A-2.31) 


Again  using  the  tame  teaneratura  profile  end,  following  the 
procedures  as  Tims  Period  1 we  obtain 


(A-2.32) 
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APPENDIX  3 


K. 


Characteristic  Equations  for  tha  System 
Five  adjoint-eigenvectors,  W/  , W/  , V/  , W/  and  W/  of 

1 2 I *t  5 

They  arc  solutions  of 


(MT-A,l)VW'  *0 
(M’-AjjW/i  ’ 0 
(MT-AiIjW/,- 0 ' 

(MT-A,i)V//4  * 0 
(M’-AsljW/j  - o 


? 


is  ths  transpose  of  the  owtrlx  M: 


O 


o 

"Ms 

0 

0 

0 

0 

-lies 

0 

0 

0 

0 

a-fi) 

-§ 

0 

° f3as 

-Up 

Tbs  five  eigenvalues  X ^ , 
(3.6*9) 


X sad  X 
• s 


srs  glass  bj 


The  first  eigenvector  fro*  (A-3.1)  is 


matrix 


(A-3.1) 


(A-3.2) 


Eq-  (3.6.5)  - 


<*■&.-¥  0 0 

' 

Wu 

> 

0 

0 Cg  0 0 o 

Wit 

0 

- p~  o o o 

Wj 

* 

0 

0 0 •u^i>i,c3  (1-0) 

wu 

0 

-$  6 0 thi 

Ws 

0 

v J 

> 

Eq.  (A-3.3)  lop  11m  W - 0 (A-3.3) 

12 


c3w;,-^^v<3*o 


-$rW,,+  CgWis  - 0 


(A-3.3) 


••ttlag 

w„-l 


WJ  luv» 


(A-3.6) 


w«"3^ 


*>-  r 


(1-3.7) 
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(A-3.8) 


(•Up  4 us  + Cg  ) W**  ( l~0)Wis  m W,s 


and 


f'Uf+^  + Cjj'^s  - 0-W,,  (A-3.9) 

combining  (A-3.6)  - (A-3.9),  v«  obtain 

VA/.  O-0)ta  ___(,~0)PS 

vv*  m-Up^i-Cj)  (-OrUs+cjr]  * 

/ C*  .1  (A-3.10) 

( ^flRTjAr-up+t^f  7 

^-w$S%Xln- 

[ (59RT7>l»(-UetU)*C,l'1]  «-3-U) 


By  the  tan  toYtm,  the  other  four  eigenvector*  are  determined: 


1«U 


(A-3.13) 


The  component*  of 


are 


Wst 


± 

r 


Vfa-(t0fcu*+tb)(ir-i)(JRT+Pk) 

tC*-(^/H-ii9)*J02rjRT 


Tho  components  of  W/  are 

% 

W4I  * 

1 

W4S-~^- 


(A-3.14) 


The  coeponants  of  tf/^  axe 

Ws>— Wit  -W»-0 

(A-3.15)  j 

v%»“Ti^r 

| 

Since  W/  , tf/  , tf/  , W/  and  V/  are  linearly  Independent  vac  to  re,  j 

1 » 1 n < | 

we  nay  Multiply  the  vector  differential  equation  Bq,  (2.6.1)  by  U/  , j 

i i 

i 

V/  » 5?/  1 W/  end  V/  aae  oorau  rive  equivalent  aquations  In  a new  cat  j 

t ’ 1 -■.*.•  « 

of  dependent  variables.  For  1*1,  2,  3.  4,  5. 
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i at 


f ' 

“3 

' X,  ' 
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T 

It 
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P 

Is 

<4> 

us 

u 

0 

J 

. 

Is 

4 / 

(A-3.16) 


Tht  Jahoaoganaoua  Ctru  I , I , I , I and  I art  dafiaad  by  Eq. 

12  14  t 


(3.6.2).  According  to  (A-3.1),  we  have  the  equality 

W/TM-A-vyf 

vrsi-w 

W/,7V1-A3^t 

tftglv l-W 


(A-3.17) 


Thus  (A-3.16)  bn 


a 

^3 

’ I, 

T 

P 

T 

P 

»iW{T 

Is 

t*P 

U 

» 3»  4, 

3 

0 

V > 

Is 

i 

(A-3.18) 


2 • Z , Z , Z and  Z at  tha 
i i » * i 

$ •atltfyltvg  tha  following  ewidmoaa 


funedent  of  u , X,  f,  u 

t P 
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W«  than  have  the  idantltiaa 


(A-3.19) 
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(A-3.20) 


(A-3.20)  can  b«  rnwxlttan  as 

Al  #♦»#•!  (*-3.21) 


•h.  (A'3.21)  la  tba  charactariatle  rapxasantatloc  of  Eq.(3.£.l). 

this  has  baan  obtalnad  by  changing  tba  dapandant  aariablaa  og,  T,  ?, 

a.  and  4 to  «aw  aariablaa,  Z , Z , 2 , Z and  Z da  tba  transformation 
P x * i * s 

(A- 3. 16)  and  (A- 3. 19) . Tbass  characteristic  aquations  can  ba  raduead 
further  by  rirtaa  of  tba  dafiaad  charactariatic  directions  in  Kq.  (3.7.1) 
and  (3.7.3) 


IAS 


Silica  Z , i*  , Z , Z and  Z are  not  perfect  differentials , cba 
i i i «*  s 

characteristic  aquation  (A-3.22)  can  be  represented  by  the  differentials 

of  u , T,  P,  u and  <j>  along  tha  charactariatic  lines . From  the 

t P 

definitions  in  Bq.  (A-3.3)  to  (A-3.11)  and  (A-3.19)  wa  have 

(A-3.23) 


Sinn*  2 - X I,  ?B  u^,  *),  by  tha  chain  rule  ve  have 

1 • 8 P 

^ (A-3.24) 
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Similarly,  v«  have 


d£«_l^£,u9+l ¥ar*Wap  u*3-2j) 

di*m$ugda9  »$F,<JT  +WdP  *f&AUf+W**  <A*3-2t) 


d£-f&d“5  ^dT  *$**+£{**  (A-3.27) 

dZ*‘*Zda3  t|f,dT  rsf>el'5  1'|ui5i:!uC W-3.2W 


Mtitltiitinj  £q.  (A— 3.24)  - (A-3.23)  into  (1-3,12)  in  fat 
wM'*s‘h+w*  (dT),  t^ap^.Muiea-^mtiop),  (*-3.29) 

-Lw/:- 1 J« 

Wm(dUf)MrWaldT)3  (A-3.3M 

-M'-J]* 

Wt  W«S  Jb  ♦ K4»Wtb  t Mfe  (dPi  ♦ W{,  (dupVMfcPAb  UrS. ») 

-faf-IJ* 

VVU(^4»^(dTVW»WP^+'«-to«pir'Wi»W#l  (A'3-32) 

— (Wf  I J*t 

AUduj)r  tV^«D»+V\W(<«P)»-»  ^(dupVfVaWPi  (*-3.33) 
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Whs  re  the  v*ctor  product*  of  W/  , W/  , W/  , W/  and  W/  with  I 

12  1%  5 

in  given  as 

W/,T'im  W,t  I,  ♦ Wi  la  * W,|  I j + Wia  I*  * Wit  Is  • 

V$T- 1 -\A4i  L +IV«  L+WoI,+ W14 Ia-H/Wi  I, - Ki 

Ii  f VVlt  Il+Wjals  *l/V»l4  + WjsIs“K£i  (A-3.34) 
W/vI“V\4<  Ii  +W*t  It  +W»»  Is  +YVUI4.  *!A£*5l* 

I,  +W4*  I*  Mis  I3  -»Ws»  U I» 

Subs ti tut ins  Sq.  (A-3.34)  ls.ee  (A-2,24)  - (A~3. 33}  and  using  th* 

daf ini t ions  of  I , I , I , I and  I , va  obtained  th*  characteristic 

1 2 j % s 

aquations  (3.6.11)  - (3.6.17). 

Besides  the  above  five  characteristic  aquations , w*  still  require 
the  consideration  of  another,  yet  to  be  determined , usual?  the  particle- 
path  line  In  tha  particle  phaaa. 

Fran  tha  haat  aquation  of  tha  particle  phase  ve  obtained  Eq.  (3.3.9) 
rewritten  as  follows. 


IfrtT  -MMfcjffl. 


(A-3. 35) 
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for  convenience  sake  we  define  the  following  par an* tar a,  tha  in- 
homogeneous  tarn  in  Eq.  (A-3.35)  la 


m or**  o<  ^ i*  I— d w b.  -eh  ~U« 

than  M , *o  v*  dacideu  to  dreg  tha  9 ten. 
i i 

W*  now  haw*  a hatter  fan  of  Iq.  CA-3.35) , 

1 *N,  (£-*“,>  gL)  («.j» 


Again  dwfla* 


Tharafora,  wa  gat 

(4Hr  1 +N '(■&). 


(A-3.41) 


Tha  sixth  char me tar Is tic  aquation  can  ba  axprassad  as 


for  Tina  Parlod  1. 

A similar  aquation  darlvad  for  Tina  Parlod  2 Is, 

d Tw)a  + K * (dT Kt> 


(A-3^2) 


(A-3.43) 


In  summary,  tha  final  foxns  of  tha  charactaristic  aquations  darlvad 
mb ova  ara  shovn  In  tha  following. 

Tha  right- running  charactaristic  aquation  In  tha  gas  phasa 


Ki,  (<*«*)»  ♦ Kic dp\  tKjpCdopHK^),  -fOt 


U-3M) 


Tha  eoafflelanes  of  tha  abova  aquation  ara, 


Ki-i 


(A-3A5  ) 


( A"3  46  ) 


(A-3.47) 


130 


Th*  l*ft-rusning  ch*r*cterietlc  *qu*clon  is  th*  |u  ph*a* » 

Kug(dtt3)a+Kf  (dPb, + K^(dUj)*+Kj (d0i*K£,£t  u-3  .50 ) 

The  coefficients  of  th*  shov*  *qu*tit«  m( 
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(A-3.51) 


Ki-i 

Kp  ‘ikj 


(A-3.52) 


f(-Up-^,-c>»?Tj.ptJgr  c»  ,ii 

ICcT. ,- a,»a5 - c}J‘Ja OlPvfrU^ u,-c/ 'Jf 


(A-3.54) 


(T-Tn^+fr  fb(hch*m 

~C*T  + ,—gP 

*'  * mt  JTlWu&iPcfr  PLc'-tur^s-Csfjp 


t 


j ysRUf-uptctj  qj* 1 

Hi 


J;  ( mcv-upra,-a>j 


(A-3.55) 


Th«  $M«ou«*^cth  characteristic  •qurtlon  1*  tu*  fas  phaaa 
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(A-3.56) 


Th*  coefficient*  of  the  above  aquation  are. 


U-3.57) 


(A-3.58) 


f'u*  [C,-C-Up»  Uj)*J«S»  JRT 

K*  FTu rfi&JtfrW 

IJa^-  f>r(T-T»hi>/;[h 


(A-3.59) 


U-3.60) 


-tx* “fifty'  ^(M-0-crf^> 

y(At  -hr(T-T*)+fttfhQ*n>-CJ  '♦ 

4A0Dr(H)Cf^-iW  / A»n>c§CW  i 

^(-^u^5s»r  1 im m&fiww 


(A-3.61) 
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Tha  right-running  charactariatic  equation  In  the  particle  phase  la 


Kup(dttf^a+ 

(A-3.62) 

The  coefficients  of  the  above  equations  are, 

Kup“i 

(A-3.63) 

k;  ~t&) 

(A-3.64) 

!/■*_  A*9Dt  AsKc 

&(i-0)  a-#) 

(A-3.65) 

The  left-naming  characteristic  equation  in 

the  particle  phase  is 

At 

(A-3.66) 

The  coefficients  of  the  above  equation  are, 

K^-  1 

(A-3.67) 

K*“  ) 

ki- 

0-3.69) 

The  panicle-path  characteristic  aquation  in  the  particle  phase  in 
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either  (Tlaa  NtRUMl  1), 


U-3.70) 


(dWk^TJj.K 


n 

Wr. 


At 


or  (Tim  Period  2) , 

K{J.t(dTw)r  » K*  fdT)e 


(A-3, 71) 


lb*  coex'ricient*  of  ehfc  above  equation*  ere. 
For  Tim  Period  1: 


K?„,  ‘ i 

(Ar3.72> 

^Ti  U^ferf.-JMhrirry.J 

(A-3. 73) 

i,«  f kr<Tf,-r.btf.hrn-TnSj 

rN,“.“  l / 

(A-3. 74) 

• 

Ti»e  Pexiod  2: 

K^'! 

(A-?. 75) 

Kr* -(,§%•) 

(A-3.76) 

(A-3. 77) 

. m 


APPENDIX  4 


Determination  of  Eigenvalues  for  the  Governlni  Equation* 


The  purpose  of  this  appendix  la  to  Indicate  the  algebraic  steps 
to  obtain  the  eigenvalues  given  by  Eqs.  (3.6.4).  The  governing 
equations  are  t reasoned  into  a matrix  form  given  by  Eq.  (3.6.1). 

In  order  to  find  the  eigenvalues  of  the  system  wt  have  to  set 
| H - X I | • 0,  where  H is  given  by  Eq.  (3.6.3)  and  1^  is  the  identity 
matrix. 

Nov  we  have  the  following  matrix. 


0 


0 


0 0 


(A-4.1) 


0 


CE 


0 -(uptA) 


0 (b0) 


Be  them  perform  the  following  alaosntary  rev  end  coIuk  rperatlon 
to  arrive  to  c better  form  to  eperat*. 


Interchange  Colusa  1 with  2 and  expand  we  get 


-(■O^+A.) 

^ 0 

& 

z£jte 

8 

'(U3  ^ ' p?RT 

0 

0 

0 

a* 

U-*) 

0 

0 </-0) 

-(Up*  A) 

Expand  Column  1 we  have 

~(«4+K) 

-O-0)P* 

a 

JHpRJ 

0 

-(■Up*  A.) 

- * 

0 

0-0; 

0 

1-0 

'-a 

0 

-gP 

0> 

0 

- <«*♦*) 

c4 

0-0) 

o 

a-t) 

(A-4.2) 


(A-4. 3.1) 


(A-4.3.2) 


130 


from  (A-4.3.1)  va  havt 


-,UfrK) 

-(<*4+K)* 

(l~0)  -Up**-) 

» ) (A-4.4.1) 

fro*  (A-4.3.2)  «•  have 

- ifT) 

nuj  *ay{ 

(\-p)  -(«p+x) 

i 

Expand  Eq.  (A-4.4.1)  and  (A-4.4.2)  wa  gat 

c *{(uPf  aAcaJ*o 

(A*4«5) 

% 

♦A)V«4y*A)CStJ[UipfA)ifi*>«>  (A-4.6) 

J 

* 

;i 

The  zoota  of  Sq.  (A-4.6)  art  1 

* 

I 

A,  + 

Aim~  149-Cg) 

A4--  <*3 
Aim~  (<£p+c) 

A«»- <Cip~C) 

fzota  tha  parftieia  haac  aquation  wa  obtain 

A * — - Up 


159 


(A-4.7) 
(A-4. f) 
(Ar4.9) 
(A-4.10) 
CA-A.ll) 


vfa«n 


(A- 4. 13) 


c,  »p««d  of  sound  In  particls  phaa«,  an  givon  by  Eqs.  (A-6.1)  or 
(Ar-6.2) . 


JE1E£2_* 

Finlto-Diffaroaca  Ecuationa  for  the  ft—  Mmi  Equation 
This  appendix  1*  intended  to  show  a typical  way  of  expreaalng  tba 
governing  aquation  into  ita  finite-difference  fora.  Tba  gas  aaaa 
aquation  was  chaaan  bare  tor  simplicity.  First  tba  divergence  fora 
of  tba  fovarnlag  equations  wars  obtained  froa  Bq.  (3.4.1  ).  For  clarity 
we  rewrite  it  in  tba  following 

2^+l!^).A.Pwr.  u-j.i) 

d L d* 


Expressing  tba  above  aquation  according  to  Eqt.  (4.2.4)  and  (4.2.5} , 
we  have  for  the  first  step  predictor  and  corrector  as  follows. 

Predictor: 

(4-5  >2) 

and  fcrractor: 

ta[tppU4)f,r(p0a)£l‘]}tit[(i-*XA»frri)!  (*-5.3) 

Zn  tba  corrector  calculation,  tba  spatial  darxvataa  and  tba  ln- 
bnmjiiiTTT  tarns  are  averaged  be  twain  tba  J **  and  j ♦ 1 level  tarns 
aslag  a waigbt  factor,  I,  whom  0 £ 9 <,  1.  Tba  superscript  4,  indicates 
a pradl  atod  quantity. 


Ul 


Similarly  tha  aacond  stap  pradictor  and  corractor  ara  axpraasoU 


according  Co  Zqa.  (4.2.8)  and  (4.2.9),  rasulting  in  two  analogous 
•quatiora . 

Pradictor: 

(A-5 . 4) 

*Z  4t  (A,  P,r.)f 

Corractor: 

f fi-e)(o=0u9  i„  - ippu3)U  ] 

+ flfcP0uj)^'tfttU3  tH >% -<P0Uf£ ]} 

+At  (( i-8) (Asp,. r + 0 (A*# *(AsPrrM  (a-5 .3) 

In  tha  corractor  calculation  thraa  lavala , j th , J + 1 and  j + 2 ara 
uaad  togathar  to  avaraga  tha  apatial  darlvativaa  and  lnhoao ganaoua  tarna 
In  Eq . (A-5 . 1) 

Zqa.  (A- 5. ?)  - (A- 5.5)  ara  tha  typical  saquanca  for  axpraaalag 
tha  conaarvatlon  aquatlona  ualng  our  nuxarical  achama.  Othav  con- 
aarvatlon  aquatlona  ara  codad  In  a similar  aanaar. 


m 


APPENDIX  6 


The  Ralationahipa  of  Intergranular  Street  and  the 
Sp — J o*  Sound  In  Granular  Md» 

This  appendix  la  aainly  to  derive  a t (tha  intragrfcnular  street) , 

P 

expression  in  tana  of  c (the  speed  of  sound). 

We  proposed  tha  epaad  of  aound  in  tha  solid  propellant  aa  a 
function  of  porosity  and  la  described  by,  (Tig.  26) 


0 £ 0c 

Crafe  »♦  0>0« 


(A-6.1) 

(A-6.2) 


ia  obtained  from  Soper's  (1C, 98)  work  and  it  is  a dimuaionlaaa 
exponent  for  tha  dacay  of  c after  +c  ia  reached. 

The  velocity  of  propagation  ia  given  in  Inference  (59)  aa  • 


(A-6.3) 


Mem  lntagrata  Xa.  (A-6.3)  for  tha  packed  bed  region. 


to*  6 £ 

After  integration  and  rearranging  of  tarn  va  gat 


(!"♦)  ~ ?(!“♦_)*  !!(♦  ■♦)  m Speed  of  Sound 


“ g [o  - 0.  HP)  t ( I - 1>)  T>]' -p,  c *0 l( 


(A-«.5> 


For  packad  bad*  T la  axpraaaad  as  follow*, 

P 

-xrmJfc$>9§i (A-0)--fcf-P  (A-‘-w 

for  $ <*c. 

Is  th«  caaa  of  dicparsad  parti  daa , w»  again  integrate  Sq.  (A-6.3) 
fro*  A to  1 ve  haws, 

gf(i-0Kr.)^-P,c^'  t*'*  (a-6-7) 

Stow  define  r,  ! • +e  and  subaticsta  is  Kq.  (A-b.7) 

9 0"f  *»->  mP?Crvt(  “JJ^C  (A-A . 3) 


Af tor  lew  algebraic  —ipwlatiwa  va  gat  tha  t expraaslos  for 

P 

a fluid! aad  bad. 


r - -Jgc1  f ) ^tHVl 
> U)~W£  u * -* 

If  A lx  app  Turning  1.  tha  bracks  cad  tarn  of  Sq.  (A-A. 
approaching  aero,  aa  a result  va  have 


(A-4.9) 
»)  i* 


When  $ Is  largar  than  $c  by  s small  amount,  c2  goss  to  zero 

vary  rapidly  sines  it's  an  exponential  dacay  function.  Tharafora, 

T Is  almost  zero . 

P 

As  a rasult  It  Is  a good  assumption  to  have  tha  following  aquation. 


O-0jTp-O  for0>A 


(A-6.11) 


Whara  Fig.  27  shows  tha  relationships  batwaan  tha  intergranular  strata 
and  porosity,  In  granular  bads,  and  H Is  tha  Haavlslde  function. 


APPEOTIX  7 


,t)  is  positive  (+) ; 

up(xL  ,t)  is  positive  <-*■); 

Kon-fluidi*«d  end  no  burning 
st  left  boundary; 

Characteristics  considered  : II,  V 


% 


✓ v 

/ £ 
> . VI 


LBC2  : 

u (x,  .t)  is  positive  (+ ) ; 

E ** 

up(«L,t)  is  negative  (-){ 

Fl4idls«d  at  left  boundary; 

Characteristics  considered  : 
II,  IV  or  V or  77 


IV  or  V 
or  VI  / 


LBC3  : 


® (a^.t)  is  negative  (-); 
up(*L,t)  is  positive  (♦); 

Fluidised  at  left  boundary; 
Character is ties  considered  t II,  I 


LBC4  : 

°g(XL  ,C)  la  nugativa  («); 

Up<XL  ,t)  la  aagativa  (-); 

Iluidizad  at  left  boundary; 

Characcarlaeica  conaidare.d  : 
II,  II’.,  VI  or  V or  VI 


LBC5  : 

u^(s^,t.)  ia  nagaelva  («) : 

V*L  „t)  it?  poeitiva  (+) ; 

Non-fluidizad  at  loft  boundary; 

Charact aria t lea  couaidac&d  t 
II,  III,  V 


lBC6  $ 

u^Cx^t)  iB  ********  C— 3 1 

^(x^t)  la  uagatlva  v-)i  j 

Nour-fluidisad  at  left't  boundary ; ( 

Char ac car iar lea  eonaldarad  : ] 

XI,  III,  V,  VI  i 


LBC7  : 

up(*L»t)  1 m a«|>civ«  (-); 

”on-£luidlz«d  at  left  boundary; 

Characteristic*  considered  : II, 
V,  VI 


LBC8  : 

u^Cx^.t)  is  positive  (+) { 
u^O^.t)  i»  positive  (+); 
Fluidised  at  left  boundary; 
Characteristics  considered  : II 


This  appendix  is  to  chov  the  various  cases  m considered  for  the  left 
boundary  conditions  by  varying  the  floe  directions  of  u^Cx^.t)  and  ^(x^.O 
sad  tbs  fluidisation  eruditions  at  the  loft  boundary*  These  eight  sets  of 
left  boundary  conditions  hove  been  inplsnnntsd  la  tbs  conputer  progrsn 
individually* 


U» 


NOMENCLATURE 


Lttia^SgbgjjL 

A 

A 


A 


k' 

1 


b 


Satcrineioa 

Pr*-*acpon*ntial  factor  for  burning  rat*  law 

Crasa-aactional  araa  of  th*  cylindrical 
coKbusti.cn  cbaabar 

Specific  aurfac*  araa  of  th*  granular  pro- 
p*li«atat  tha  surface  exposed  to  fluid 
par  unit  volvaa 

Tb*  intersections  of  eh*  1,  II , III,  IV , V, 
VI  charactAriatics  with  tb*  x-axis  aear  th* 
laft  boundary  of  tha  granular  b«d 

Tb*  an  tar  Motions  of  tb*  I,  II,  XII,  IV,  V, 
^I  character Aaticc  with  the  x-axls  near  th* 
right  boundary  of  tb*  granular  bad 

Co-velum*  of  the  Claus  i-us  aquation  of  state 


B 


Bovndery  point  at  tha  laft  and  of  tb* 
granular  bod 


B' 


e 


Boundary  point  at  tha  right  *e&  cf  the 
granular  bad 

Syaad  of  sound  in  panic]*  phaa* 


C 


d 


D*tarr*nt  concentration 


t,  Co*£uetloa  eh amber  length 

Spaclfle  boat  at  eonatant  volume 


Swtl Lie.  L«at  at  cerastaat  pr^'err* 


it 

dt 


Slops*  of  th*  that  act*  viatic  <kj  cat  ictus 


Tb*  total  <]£k|  fore*  between  tha  gas  and 
partial*  pisaoos  *»  tb*  antrsneo  region 

Brag  force  duo  to  porosity  pradiant 


D Tb*  total  dvag  force  batwean  eh*  gao  and 

$artlcl«  rb****,  • 3„  + d 

▼ p 

5 drag  lores  acting  on  gaeat  by  the  particles 

per  ante  w*tt*d  era*  of  particles 


tagrmgian  tiat  derivative  wh*a  obaarvir 
follow*  th*  notion  cf  a part  id*,  9 . I 

It  ' 


Met* r Symbols 
E 


h . 
ch%o 

hign 


K 

I 


I 


m 


I 

,a 

J 

k 

k 

P 

X 


. K 


ia 

mm 

tt" 


Inscription 

Total  scored  energy  (Internal  plus  kinetic 
energy)  par  unit  mass 

Gravitational  constant 

Tha  avorafs  couvsction  heat-transfer  co- 
at' fi  clan  f eves.'  palls t a 

Enthalpy  of  eha  propallant  gas  at  flans 
temporal  ire 

Enthalpy  of  the  ignitar  gaa 

The  svnrugs  radiation  hast  transfar  co- 
efficient ovar  pellets 

Tha  total  haat  transfar  coefficient,  the  ;ait- 
aurface  conductance 

The  rata  of  change  of  unit-surf aca  conductance 
in  foil  O' ring  a civen  particle 

Baaviaidc  function 

Inhomogeneous  vector  for  tha  governing 
attentions 

Identity  matrix 

Component t of  tha  inhomogeneous  vector  I 
Haat  tj  wtrk  conversion  factor 
Tharaal  conductivity  of  gases 
Thermal  conductivity  of  particles 
Coefficient  of  charactsrlatic  aquations 
Eromlva  burniag  constant 
Mass  flow  rata  of  Ignitar 

Mass  flew  rata  of  gaa  into  tha  granular  bad 

varitble  matrix  defined  ia  Eq.  (3.6.3) 

Prmmaura  aapomant  for  tba  burniag  rata  law 

Buafear  density  of  apharical  pa lists  in  tha 
granular  bad 


172 


Pressure 


Prendtl  number 

Wetted  perimeter  between  the  chamber 
wall  and  gas  phaae 

Wetted  perimeter  between  the  chamber 
wall  and  particle  phase 

Kata  of  conduction  haat  transfer  par  unit  area 


Bata  of  haat  loss  to  tha  chamber  wall  par 
unit  spatial  volute  of  gas-particle  system 

Badial  distance  from  tha  center  of  a spherical 
particle 

Burning  rata  of  the  solid  propellant 
Badius  of  pellets 
Gas  constant 

2r  pd|u  - uj 
Baynolds  number,  1 "*  1 1 •- 

Time 

Temperature 

Ablation  temperature  of  the  solid  propellent 
Adiabatic  flan*  teaperature  of  pellets 
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